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ABSTRACT 


The effect of tidal friction on the inclination of the lunar orbit to the earth’s 
equator for earth-moon distances of less than 10 earth radii is examined. The 

results obtaineu bear on a conclusion drawn by Gerstenkorn and others which 
has been raised as a fatal objection to the fission hypothesis of lunar origin, 
namely, that the present nonzero inclination of the moon's orbit to the ecliptic 
implies a steep Inclination of the moon's orbit to the earth’s equatorial plane in 
the early history of the earth -moon system. This conclusion is shown to be valid 
only for particular rheological models of the earth. In the case of a viscous 
earth , the results indicate that the problem of wrenching the moon out of an 
equatorial orbit into an inclined orbit to account for the present tilt of the lunar 
orbit to the ecliptic must be faced in the accretion theory of the moon’s origin 
and possibly the capture theory, as well as in the fission theory. In this respect 
all three theories are on the same footing. A solution to the inclination problem 
is presented. 

The treatment of tidal friction adopted here employs the approach of George 
Darwin and pursues his suggested solution to the inclination problem in great 
detail. The earth is assumed to behave like a highly viscous fluid in response to 
tides raised in it by the moon. The moon is assumed to be tideless and in a 
circular orbit about the earth. The equations of tidal friction are integrated 
numerically to give the inclination of the lunar orbit as a function of earth-moon 
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distance. It is found that if the radius of the lunar orbit is greater than 3.83 
earth radii, then the inclination of the moon's orbit to the earth's equator will 
increase if the moon is perturbed from an equatorial orbit, provided the earth's 
viscosity is greater than 10^® poises. The present inclination of the lunar orbit 
to the ecliptic can be explained if the moon's orbit is perturbed about 3° out of 
the equatorial plane at 3.83 earth radii, provided that the earth's viscosity is 
not less than 10^® poises. It is also found that if the viscosity is large (greater 
than 10^® poises), then, under certain conditions, the, radius of the moon's orbit 
may actually decrease temporarily, and then increase; and further, that an 
upper limit can be placed on the inclination of the lunar orbit to the earth’s 
equator when the moon is 3.83 earth radii distant from the earth, regardless of 
the moon's prior history- 
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PREFACE 


Readers unfamiliar with tidal friction should find Chapter I, 
Section A and Appendix A of some value. A list of important 
quantities for this work is given in Table 4. A list of correc- 
tions of misprints in Peter Goldreich's important paper "History 
of the Lunar Orbit" is given in Appendix F. Page numbers of 
the reference "Darwin (1880)" refer to Darwin's paper as it ap- 
pears in Scientific Papers by George Howard Darwin, Volume II, 
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INTRODUCTION 


Men have speculated about the tides for centuries. An ancient Chinese 
scholar suggested that the earth lived, that the ocean was its blood, and that the 
tides were the beating of the earth’s pulse (Darwin 1962, pg. 76). An Arabian 
scholar explained the rising of the tide as being caused by the heating of the 
ocean by sunlight and moonlight (Darwin 1962, pp. 77-79). Others also suggested 
that the tides were somehow caused by the sun and moon (Darwin 1962, pp. 79- 
85), but it remained for Isaac Newton to advance the correct explanation for the 
cause of the tides (Newton 1966). Newton realized that the lunar tides were 
caused by a combination of gravitational pull and centrifugal effects which would 
make the water in the oceans collect on the sides of the earth directly under and 
directly away from the moon, thus giving the earth a bulge. A similar argument 
holds for the solar tides. 

Newton's theory of the tides was carried forward, notably by Bernoulli, 
Laplace, Darwin, and Kelvin (Darwin 1962, pp. 86-88) to explain the rise and fall 
of the oceans on the earth. Their efforts culminated in the work of Doodson and 
Proudman (Doodson 1958). 

George Howard Darwin, son of the famous Charles Darwin, considered not 
only the problem of the tides raised on the earth by the moon, but also a more 
subtle problem: the action of the tides on the motion of the moon (Darwin 1880). 
He included in his investigations not only ocean tides , but also tides raised in the 
bulk of the earth as well; these latter tides are called earth or body tides. 
Darwin recognized that friction attending the tides, whether they are raised in 
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the oceans or in the earth, would have profound effects on the moon's orbit. In 
fact, tidal friction dominates the secular chaise in the moon's orbital elements. 

Darwin assumed the earth to be a homogeneous, incompressible viscous 
fluid in which tides Were raised by the moon; the moon Itself was assumed to be 
a point-mass. He expanded the tidal disturbing function in a Fourier series and 
integrated the equations for the secular change of the moon's orbital elements 
backwards in time in an attempt to uncover the past history of the moon. Darwin 
found that the moon orbited very close to the earth at some time in the distant 
past. He speculated that the earth and moon were once a single primitive body, 
and that resonance vibrations set up in the body by the sun caused the body to 
fission, thus throwing the moon into orbit about the earth. Tidal friction then 
caused the moon to move away from the earth to its present distance. 

Jeffreys (1930) found that dissipation in the primitive body would be so 
great that the vibrations would be damped out, making it impossible for the 
moon to be torn out by the action of the sun. The fission theory of the origin of 
the moon fell out of favor. It has been reproposed in more recent times by 
Cameron (1963), Wise (1963), and O'Keefe (1969). 

Modern interest in tidal friction was rekindled by Gerstenkorn (Alfven 1963), 
who invoked tidal friction in a new hypothesis of lunar origin: capture of the 
moon. Gerstenkom's analysis lead him to propose that the moon was once an 
independent planet in an orbit which carried the moon close to the earth. The 
tidal interaction between the moon and earth captured the moon in a retrograde 
orbit, which subsequently flipped over into the prograde orbit we see today. 
MacDonald (1964) si 4 >ported a many-moon hypothesis of lunar origin to 
overcome a time-scale difficulty in tidal evolution. Singer (1968) inves- 
tigated the problem of prograde capture, while the analysis of Goldreich (1966) 
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lead him to favor accretion of the moon from a swarm of particles in orbit about 
the earth. 

The problems associated with the intimately connected questions of the 
moon's origin and its orbital evolution are seen to be of surpassii^ interest 
today. We will investigate here one of the problems of the early history of the 
moon; the inclination of the lunar orbit. 



CHAPTER I 


TIDAL FRICTION AND THE INCLINATION PROBLEM 

A. Qualitative Aspects of Tidal Friction 

Some of the qualitative aspects of tidal friction will now be examined; for 
fuller discussions see Goldreich (1972); MacDonald (1964); and Jeffreys (1962). 
We will begin by dealing with a simplified picture of the earth-moon system. 

The earth and moon are assumed to be the only two bodies in existence, with the 
moon orbiting the earth in a circular orbit Ijdng in the plane of the earth's 
equator. In addition, the moon is assumed to be perfectly spherical, and the 
earth to be without atmosphere or oceans, so that we are concerned only with 
body tides in the earth. 

Figure 1(a) shows the case where the earth exhibits no internal friction. In 
this case if the earth behaved like a solid it would be perfectly elastic; if the 
earth behaved like a liquid, it would have no viscosity. The tidal forces acting 
on the earth cause it to bulge along the line joining the centers of the earth and 
moon. The part of the bulge nearest the moon is raised by the pull of the moon's 
gravity, which is greatest on the side closest to the moon. The part of the bulge 
opposite the moon may be thought of as being thrown out by the centrifugal force 
associated with the motion of the two bodies about their common center of mass. 
In this case, there would be no evolution of the moon’s orbit. The moon would 
still revolve about the earth in a circular orbit, with only a slight change in the 
earth's gravitational force from its value for an undistorted earth. 

The situation changes, however, when friction is present in the earth; this 
case is shown in Figure 1(b). In the simplest picture, the action of tidal friction 
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makes the axis of the bulge swing away from the line joining the centers of the 
earth and moon and reduces the size of the bulge. From the viewpoint of inertial 
space (a frame fixed with respect to the distant galaxies) , the bulge may be 
thought of as being carried around by the earth's rotation. Note that in this case, 
where the angular velocity of the earth is greater than the angular velocity of 
the moon (relative to inertial space), the bulge leads the moon. The behavior of 
the bulge may also be understood from the viewpoint of an observer standing on 
the earth's equator. The observer would see the moon rise in the east and set 
in the west because the angular velocity of the moon relative to the observer is 
clockwise. In the frictionless case, a high tide would occur when the moon 
reaches the observers zenith. If friction is present, however, the high tide does 
not occur until after the moon has passed the zenith, since friction causes a de- 
lay. Hence to the observer standing on the earth, the tidal bulge lags behind the 
moon. If two lines are now drawn, one along the axis of the bulge and one joining 
the centers of the earth and moon, we get exactly the case shown in Figure 1(b) . 
The tidal lag angle is the angle between the two lines . 

If the angular velocity of the moon were greater than the angular velocity of 
the earth, as in the case of Phobos orbiting Mars, or if the moon revolved in a 
sense opposite to that of the earth, as in the case of Triton orbitii^ Neptune, 
then the bulge would lag behind the moon (as viewed from inertial space). 

We return to the simple system shown in Figure 1(b). The moon's gravity 
pulls on the nearer part of the bulge with greater force than it pulls on the 
farther part of the bulge, producing a net torque on the earth. This torque acts 
in a sense opposite to the earth's rotation; hence the earth slows down. By re- 
action, the bulge will exert a torque on the moon, causing the moon to "speed up" 
and move away from the earth. Thus the moon was closer to the earth in earlier 
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times. The total angular momentum of the system is conserved in this process, 
but the total mechanical energy decreases as friction dissipates the energy into 
heat. 

Unfortunately, things are not as simple in reality as shown in Figure 1(b). 
For one thing, the moon's orbit does not lie in the earth's equatorial plane; nor 
is it circular but elliptical, which means the distance between the earth and 
moon is continually changing. Also, body tides are present in the moon, com- 
plicating the tidal interaction (MacDonald 1964). Further, the sun also raises 
tides on the earth ; lunar and solar gravity act on both the lunar and solar bulges. 
The presence of the sun also creates a three body problem. The earth and moon 
are not spherical even in the absence of tidal forces; the earth is flattened by 
rotation, for example. Also, the actual shape of the tidal bulge is not necessarily 

as simple as shown in Figure 1(b). The shape depends upon the model chosen 
for the earth's properties. In general, the tidal forces distort the earth into a 

figure resembling a triaxial ellipsoid. 

The present-day earth has ocean tides and atmospheric tides as well as body 
tides. The varying depths of the oceans, the flow of tidal currents, and the 
irregular shape of coast lines make the ocean tides quite complex. The ocean 
tides may be responsible for most of the dissipation of energy (see below). The 
atmospheric tide is an observed semidiurnal variation in atmospheric pressure 
caused by solar heating and not solar or lunar gravitational forces. This tide 
lags behind the sun as viewed from inertial space, so that the gravitational 
torque on the atmospheric tide tends to speed up the earth. This torque may be 
comparable in magnitude to the solar ocean torque, tending to cancel it (Jeffreys 
1962). 
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Observational evidence for tidal friction conies from a variety of sources. 
Body tides are observed with sensitive gravimeters (Tomaschek 1957) from 
which the lag angle may be deduced (MacDonald 1964). The tidal bulge can be 
observed by its perturbing effects on the orbits of earth satellites (Newton 1968). 
Celestial observations, both modern and ancient (Newton 1969), reveal the secu- 
lar acceleration of the moon and the deceleration of the earth's spin. Agreement 
between the different methods is rough, but they indicate the following (Goldreich 
1972): the present-day lag angle in the simple picture of Figure 1(b) is between 
2° and 3®, with energy being dissipated at the rate of ~2.6 x 10^® ergs/sec. The 
moon is moving away from the earth at the rate of 3 cm per year, with the 
earth's daily rotation period slowing down by 2 x 10"^ seconds per year. The 
work of Miller (1966) suggests that two thirds of the energy dissipation takes 
place in the shallow seas, but the figure is uncertain and the actual seat of most 
of the dissipation (whether in the oceans or in the earth) is unknown. 

Remarkable evidence for tidal friction in the distant past exists in the form 
of daily growth bands found in fossil coral and shellfish. The work of Wells 
(1963) on fossil coral suggests that the year was about 400 days long 380 million 
years ago, which is consistent with the current rate of slowdown of the earth's 
rotation. A constant rate of slowdown over this time period has been called into 
question by the fossil evidence found by Pannella et al. (1968), however. 

It should be mentioned that tidal friction is important not only in planet- 
satellite systems but also in sun-planet and binary star systems as well. 

B. The Inclination Problem 

We will investigate possible early histories of the inclination of the lunar 
orbit. The problem of the inclination in the early stages of the moon's history 
has been cogently summarized by O'Keefe (1972). 
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Goldreich (1966) investigated the history of the lunar orbit using the 
assumptions of a circular orbit and weak tidal friction; this latter assumption 
entails either low viscosity or imperfect elasticity in the earth. Goldreich found 
that if the moon were 10 earth radii distance from the earth the inclination of 
the moon's orbit to the earth's equator would have been about 10°; at closer dis- 
tances the inclination would have been even higher. 

This result appears to rule out the fission theory of the moon's origin, 
since if the rapidly rotating primitive body fissioned, the moon would be thrown 
into an equatorial orbit around the earth, in contradiction to Goldreich's findings. 

Darwin (1880), making assumptions similar to Goldreich's, came to much 
the same conclusion; but Darwin originated the fission theory. How could 
Darwin believe in a contradiction? 

The crux of the matter is the assumptions that are made in modeling the 
properties of the earth. Goldreich and other modern investigators (MacDonald 
1964; Gerstenkorn (Alfven 1963) ) considered the effects of weak tidal friction, 
as did Darwin; but Darwin went on and examined the effects of high viscosity: 
strong tidal friction (Darwin 1880). 

Darwin found that if the moon were perturbed slightly out of an equatorial 
orbit, then, under certain conditions, the tidal forces acting on the moon would 
cause the inclination to grow. An earth which had a high viscosity in its early 
history might then solve the inclination problem. Darwin seized upon this as 
the answer and did not investigate the matter further. 

We will follow Darwin's treatment in assuming a highly viscous earth in the 
early stages of the earth-moon system and examine the inclination in more de- 
tail; specifically, under what conditions the inclination will increase for small 
initial perturbations . 
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C. Proper Planes and the Inclination Problem 

Let us now examine the inclination problem from the aspect of the proper 
planes of the moon and earth. These planes were discovered by Laplace (1966) 
and were used by Darwin (1880) in his treatment of tidal friction. 

Laplace found that the plane of a satellite's orbit about an oblate planet 
tended to maintain a constant inclination to a plane which he called the proper 
plane of the satellite. We will follow Darwin and call this inclination J. 

The proper plane lies intermediate between the plane of the planet's orbit 
about the sun (the ecliptic plane for the earth) and the invariable plane (the plane 
perpendicular to the total angular momentum vector of the planet-satellite 
system) . The angle between the ecliptic and proper planes Darwin called J/ 

If a satellite orbits far from a planet, the sun's influence dominates the inclina- 
tion and the proper plane is nearly parallel to the ecliptic (J/ % 0), so that the 
satellite has very nearly a constant inclination to the ecliptic. If a satellite 
orbits close to a planet, the oblateness of the planet dominates the inclination 
and the proper plane is nearly parallel to the invariable plane. In this case the 
satellite tends to maintain a constant inclination to the equatorial plane of the 
planet (as shown later). 

These two limiting cases are referred to by Goldreich (1966), who foimd 
that the transition between the two came at a distance which he called the criti- 
cal distance. This is the distance where the torque exerted on the satellite by 
the planet's bulge is equal in magnitude to the torque exerted on the satellite by 
the sun. The critical distance is about 10 earth radii for the earth-moon system. 

At the present time the mocai is about 60 earth radii distance from the 
earth, so that the orbital plane of the moon keeps a nearly constant tilt to the 
ecliptic. The angle between the proper plane and the ecliptic is about 8" 
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(Darwin 1880) and the angle J between the proper plane and the moon's orbital 
plane is about 5°9’. 

The earth also has a proper plane; the earth's equatorial plane tends to 
maintain a constant inclination to its proper plane. The angle between the two 
is called 1/ by Darwin, with the angle between the earth's proper plane and the 
ecliptic being I. At present ly is about 9" (Darwin 1880) and I is 23®27'. 

While the orbital plane of a satellite tends to maintain a constant angle J to 
its proper plane, the orbital plane precesses in space, so that the vector normal 
to the orbital plane sweeps out a cone around the vector normal to the proper 
plane. This is diagrammed in Figure 2(a). Likewise, the earth's axis sweeps 
out a cone around the vector normal to its proper plane. Both precessions have 
the same speed and direction. 

At small distances between the planet and satellite (i.e. when solar influence 
can be neglected) , the poles of the two proper planes merge with the pole of the 
invariable plane, so that the orbital plane of the satellite and the equatorial 
plane of the planet maintain a constant tilt to the invariable plane and each other, 
as shown in Figure 2(b). Hence in this case the orbital plane of the satellite has 
a constant inclination to the equatorial plane of the planet as mentioned earlier. 

If the moon somehow formed or arrived in the earth's equatorial plane at a 
distance of less than 10 earth radii, then the orbital plane, equatorial plane, and 
invariable plane would all coincide and the inclination J of the moon's orbit to 
its proper plane would be essentially zero. Subsequently, if tidal friction did 
nothing to affect J as it pushed the moon steadily away from the earth, then at 
the present time the moon would have essentially an ecliptic orbit (J and Jy % 0). 
However, the present value of J is about 5*9'. Thus, if it is assumed that the 
moon did form close to the earth in the equatorial plane, then it must be 
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explained why the moon’s orbit now has a five degree tilt to the ecliptic and not 
a virtually zero tilt- 

To put it another way, if the present five degree inclination is extrapolated 
back into the past, the plane of the moon’s orbit would be steeply inclined to the 
equatorial plane of the earth when the moon was close to the earth. Therefore, 
those theories of the origin of the moon which postulate the moon's formation in 
the earth’s equatorial plane must explain this discrepancy. 

The effect of tidal friction in Goldrelch’s (1966) formulation on J, the in- 
clination of the moon’s orbit to its proper plane, can be extracted from his 
Figure 7, which is reproduced here as Figure 3. 

Goldreich's figure shows the inclination of the moon’s orbital plane to the 
ecliptic as a function of earth -moon distance. The inclination of the moon’s 
proper plane to the ecliptic is so small at the present distance of 60 earth radii 
that it is imperceptible in the figure, so that the moon's orbital plane appears 
to keep a constant five degree inclination to the ecliptic for this distance. The 
normals to the ecliptic, proper plane, and the cone swept out by the normal to 
the lunar orbit for this case are shown in Figure 4(a). Note that the normal to 
the ecliptic is inside the cone. 

At about 30 earth radii the angle becomes large enough to be noticeable 
in Goldreich’s figure, so that the variation in angle between the ecliptic and 
orbital planes is clearly visible and the curve branches, showing the maximum 
and minimum inclination. This situation is diagrammed in Figure 4(b). Note 
that the normal to the ecliptic still lies inside the cone. Reference to this figure 
should make clear that the inclination J of the lunar orbit to its proper plane can 
be found from Goldreich’s figure by adding the maximum and minimum 
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inclinations and dividing by 2 . The inclination of the proper plane to the 
ecliptic is then found by subtracting the minimum inclination from J. 

At about 17.5 earth radii in Goldreich's figure the normal to the ecliptic 
lies in the surface of the cone, so that the minimum inclination is zero and the 
lower branch of the curve touches the horizontal axis. This is shown in 
Figure 4(c). 

Between 17.5 and 3 earth radii the normal to the ecliptic falls outside the 
cone, as shown in Figure 4(d). The inclination in this region is then found by 
drawing a curve equidistant between two branches in Goldreich's figure and 
measuring from the horizontal axis to that curve. The angle J is half the differ- 
ence between the two branches. 

Figure 20 shows J as a function of earth-moon distance as extracted from 
Goldreich's figure by the above process (dashed curve). Note that the inclina- 
tion of the moon's orbit to the proper plane increases as the distance decreases 
below 13 earth radii. Darwin's small viscosity model gives a remarkably simi- 
lar result (dotted line; see Chapter IV). This indicates that small tidal lag 
angles cannot be invoked to drive the moon out of the earth's equatorial plane; 
if it could, then J would decrease as distance decreases for small distances, 
instead of increasing. 



CHAPTER n 


DARWIN'S APPROACH TO TIDAL FRICTION 

We will now briefly outline George Darwin’s approach to the problem of 
tidal friction (Darwin 1879, 1880). Although Darwin treats the case of a planet 
attended by two tide-raising satellites (such as the earth attended by the moon 
and sun, where the latter may be treated as a satellite of the earth), we will re- 
strict the discussion in this chapter to the moon and earth as an isolated system; 
i,e, the presence of the sun will be neglected. 

The following assumptions are made by Darwin: the earth is a homogeneous, 
viscous, incompressible sphere. Body tides are raised in the earth by the moon. 
The moon is taken to be a point -mass without rotational angular momentum. The 
tide-raising potential generated by the moon is a second degree spherical har- 
monic. The tidal disturbing potential generated by the earth is expressed as a 
sum of second degree spherical harmonics. The effects of inertia are neglected 
when solving for the response of the earth to the tide-raisir^ force. 

We will further restrict the discussion to a circular orbit for the moon 
about the earth. 

Before plunging into a discussion of Darwin's treatment we will discuss the 
effects of the earth's rotational bulge on the motion of the moon. 

Goldreich (1966) has shown in an elegant manner that the rotational flat- 
tening of the earth produces no secular change in the magnitude of the orbital 
angular momentum of the earth -moon system or in the earth's rotational angular 
momentum; we denote these two quantities by and L£, respectively. 
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The orbital angular momentum of the system about the center of mass 
of the system is 


Lj, = 

where 

M = mass of the earth 
m = mass of the moon 

0 = angular velocity of the earth and moon about the center of mass 
= distance of the earth from the center of mass 
dj = distance of the moon from the center of mass 
Now by Kepler's third law 

ya (M + m) 

j-3/2 


where 

G = universal gravitational constant 
r = earth-moon distance 

Also 





We may now write 
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Lg is easily seen to be 


Lg C n 


where 

n = rotational angular velocity of the earth 
C = polar moment of inertia of the earth. 

We may show the constancy of other important quantities by usii^ Gold- 
reich's result of the constancy of and Lg. 

Clearly there is no secular change in the earth-moon distance r if there is 
no secular change in L^. Likewise, there is no secular change in the rotational 
angular velocity of the earth n if there is no secular change in Lg. Therefore r 
and n are constant in the secular sense for the case of the rotational bulge. 

We refer now to Figure 5, which shows the angular momentum triangle for 
the earth-moon system. L^ is the total angular momentum of the system and is 
constant in both magnitude and direction because the system is isolated. It is 
clear from the diagram that if L,^, Lg, and are unchanging, then j, the angle 
between the plane of the lunar orbit and the invariable plane, and i, the angle be- 
tween the plane of the earth's equator and the invariable plane, are constant. 
Thus the rotational bulge of the earth produces no secular change in r, n, j, or i. 

The flattening of the earth does make the lunar orbit precess in space. It is 

clear from Figure 5 that the earth must precess at the same rate and in the 

same sense as the lunar orbit. By conservation of angular momentum L.j. = 

— ♦ — ♦ 

Lj, + Lg, so that the three vectors must lie in the same plane. 

It should be clear, then, that if the longitude of the moon's node N is 
measured along the invariable plane from the descending node of the intersection 
of the earth's equatorial plane and the invariable plane, N must be zero. 
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In this investigation we are chiefly concerned with secular changes in r, n, 
j, and i; hence further consideration of the effects of the rotational flattening on 
the moon's motion will be dispensed with. 

It should be mentioned that the sun also causes no secular change in r, n, j 
and i to the order of approximation carried out by Goldreich. 

Let us now investigate the response of the earth to the tide-raising force 
and the effect of the earth's response on the moon (see Appendix A for a deriva- 
tion of the tide-raising potential and the tidal disturbing function) . 

The tide-raising potential at some point (x*, y*, z*) in the earth is given by 
Equation (A-6) of Appendix A as 

3 Gm /r*\^ f , ll 

Vt = 2 —[ t ) H®-3j 


and is the first equation of §4 of Darwin (1880) with some notational changes, 
r* is the distance from the center of the earth to (x* , y* , z*) and the angle 0 
(which Darwin calls PM) is shown in Figure 21. 

If the earth were a frictionless fluid the tide-raising force would raise a 
tide on the earth, with the height of the tide o-^ being given by (Darwin 1879, 
Equation 13): 


Cr 


t 


15 Gm 
4 g r3 



where a is the radius of the earth and g is the gravitational acceleration at the 
earth's surface. The earth would clearly bulge in this case as shown in Fig- 
ure 1(a). Note that the height of the tide is inversely proportional to the cube 
of the earth-moon distance. 
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Darwin chose axes fixed in the earth and expanded ^cos^ 0 - in terms 
of the direction cosines of both (x*, y*, z*) and the position of the moon (x, y, z). 
Letting 77 and C be the direction cosines of {x*, y+, z*) and M^, Mj, and 
the direction cosines of the moon, then cos ©=^Mj + 7 jM 2 + I M 3 and we 
may write 


15 Gtn f 

= T T 7? f 


Mj + 2 2 


^2 _ ^2 m2 - Mj2 


2 + 2 7 } CM2M3 


(II-l) 


3 ^2 + ^2 _ 2^2 M 2 +M|- 2tii\ 

+ 2#^m,m3 + 2 ^ — h — ^ 3 — ■; 


after some algebraic rearrangements. 

Mj, Mj, and depend upon i and j, the respective angles of the earth's 
equator and the plane of the moon's orbit to a fixed plane, which, in the two body 
problem, we take to be the invariable plane; n, the rotational angular velocity of 
the earth; H, the angular velocity of the moon in its orbit; N, the longitude of the 
node of the moon's orbit measured from the descending intersection of the 
earth's equatorial plane and the fixed plane along the fixed plane; and t, the time. 
For example, is given by Equation (20) of Darwin (1880) as 

Mj = P2 p2 COS (X - I N) + P2 q2 cos (X + £ - N) + p2 cos (X + £ + N) 

+ Q2 q2 cos (X“£ +N) + 2PQpq [cos (X + £) “ cos (X “ £ )] 


. Q = sin (I 

with Xq a constant, and £ = Ht + e , where e is the loi^itude of epoch. Mj and 
M 3 may be written in a similar fashion. Notice that Mj is expressed as a sum 


Here P = cos 


(l‘) 


ij , p = cos i) . q = sin , X = nt + Xq 
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of terms periodic in time, whose periods depend on linear combinations of n and 
Q; the same is true of and M^. 

The combinations of Mj, Mj, and M 3 that appear in Equation {II-l) (M^ 
m2 - M = 

, etc.) can also be written as sums of simple harmonics whose angular 

2 

speeds are linear combinations of n and (Table 1 gives the total number of 
angular speeds which arise.) For instance, after considerable work Darwin 
shows that Mj Mj may be written (Darwin 1880, Equation 25): 

“1^2 = + 2772^2 + 

(H-2) 

- jr* e“ 2 '^ (X- «) _ 2772 k 2 e-2«/^X _ ^4 ^-2/^ (X+ 


Here 77 = Pp - Qq e+’^N; /< = Qp + Pq e*^N; 77 = Pp - Qq 6 "’^^. 

K = Qp + Pq e"'^*^; and 6 = t + N. Darwin put the sines and cosines in 
exponential form for convenience in later work. 

Equation (II-l) could now be written 


cr 


t 


IS G m a2 
"4 g J.3 




+ 2772^2 62-^ X + 


■1 


by substituting in it the complicated expressions for Mj Mj, etc. 

The equation above gives the displacement of the earth's surface when the 
earth is composed of a frictionless fluid. What we now wish to find is the ex- 
pression for cTj, when friction is present inside the earth. 

It is assumed that the effects of friction are such that each simple time 
harmonic that appears in the expression for o-^. is multiplied by a factor to re- 
duce the amplitude of the harmonic, and its phase is altered by a certain lag 



angle. For example, Mj now becomes in the presence of friction (Darwin 
1880, Equation 33) 
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FT 


[2(X- ^)- 2 f,3 + F 2772 ^2 e’^ t2X- 2f] 


+ Fj K* [2(X+ 6)- 2f2] _ Fj 7t;4 [2(X- ^)- 2fj] 

- F2zr2/<2 t- 2 X + 2f] _ f^ ^4 g- [2 (X + «) - 2 ( 2 ]^ 


Fj, F, and F^ are amplitude factors and 2f^, 2f, and 2fj are the respective 
phase angles. (Table 1 gives the amplitude factors and phase angles for all the 
speeds.) 

Darwin calls the above expression becomes 3(2 - , etc., 

so that now in place of (II-l) , we have 




(D-3) 


, 3 + 7,2-2 3(^ + ^^-222 

2 3 3 


as the equation for the earth's surface in the presence of friction (Darwin 1880, 
Equation 30). 

The exact values of the amplitude factors and phase angles depend upon the 
model chosen for the earth's properties. The model we are interested in is the 
case where the earth is a fluid exhibiting Newtonian viscosity. Darwin (1879) 
found the amplitude factors and phase angles for this case, which are given by 
the following relations: 
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19 V 

tan [ lag angle] = [ angular speed] x 2go.p 

[amplitude factor] = cos [lag angle] 

where p is the density of the earth and v is its viscosity. 

An example of the above relations is 

IQ -u 

tan 2 f = 2 (n - 0) - 

1 ' 2 g a /O 

Fj - cos 2 fj 

for the angular speed 2 (n - D) . 

Now that cr^ has been found, our next step is to find the tidal disturbing 
function acting on the moon. It has been derived in Appendix A and is given 
by Equation (A-15): 

= ~7tG pa (|r) (H-4) 

where r' (which was called A in Appendix A) is the earth-moon distance, a' and 
/?' are the longitude and colatitude of the moon in the earth-fixed frame, and p is 
the earth’s density. Primes are placed on the variables for reasons discussed 
below and in Appendix A. Note that implicitly depends on r~^ (see Equation 
n-3), so that the disturbing function is proportional to the inverse sixth power of 
the earth-moon distance. 

We now wish to know how the disturbing function changes the orbital ele- 
ments of the moon. Here Darwin uses the Lagrange equations for the time 
derivatives of the osculating orbital elements. These equations are derived e.g. 
in Brouwer and Clemence (1961). Darwin uses four of the six equations in his 
1880 paper, which we reproduce in his notation (Darwin 1880, Equations 1-4): 
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dc 

20.c^ 

br 



dT 

G (M + m) 




de 

Q. c 

1 - B R Vl - /B R 

B R\‘ 


dt 

G (M + m) 

e Be e \B e ^ 

B 77 / 


11 = 

G c 

1 f 1 br 1 I 

/br 


d t 

G (M + m) 

_ ^2 j B N 2 ^ ' 

iBe 

'd 77 / 


dN 

fi c 

1 br 

dt 

G (M -(■ m) 

yi - e2 ^ J 


The only quantities not defined by us thus far are c, the semi-major axis of 
the orbit; e, the eccentricity; tt , the longitude of perigee (not to be confused with 
77 used elsewhere in Darwin); and R, any disturbing function in general. Since 
we are interested in the effects of the tides on the moon, we set equal to R. 

Darwin alters the form of the equations to make them more convenient to 
use for his purposes. For example, in the case of a circular orbit the first and 
third equations become (Darwin 1880, Equations 11 and 13): 


1 ^ 

k dt 3 e 


1 1 3 W 1 . BW 

k dt sin j B N ^ 2 ^ Be 


/ c 

where now ^ = /— j , with c^ being some reference distance, and c is now the 

c 

radius of the orbit; k is ^ ^ c^ , with being the angular velocity of the 


moon at the reference distance c^; and W is ^ R, with C = Ma 

8 


2 _ 
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The next step is to express W in terms of i, j, N, and e and evaluate the 
derivatives of W to find the time rate of change of the moon's parameters. 
Before we proceed, however, we must heed the warning given in Appendix A not 
to confuse the moon's parameters as they enter in the role of tide -raising body 
and in the role of tidally disturbed body,, even though here the bodies are one 
and the same (the moon) . Let us therefore follow Darwin and place primes on 
the parameters of the disturbed body (which we have already done in Equa- 
tion II-4). Our expression for W becomes (Darwin 1880, Equation 31); 

T r' r X '2 - Y '2 3(2 - ^2 

W = — [^2X'Y' + 2 2 2 ^ 


+ 2X’Y'3(2 + ^ 


3 X'2 + Y 


• 2 _ 


2 Z'2 3(2 + ?12 - 2 22 


where X', Y', and Z' are the direction cosines of the moon in its role as dis- 
turbed body, r = with (with a similar expression for r' )> 

j 2 g ® 

and g = -g- 

X' Y' is the same as Equation (11-2) save primed variables replacing the 

X*2 ^ yi2 

unprimed variables; and similarly for ^ , Y' Z' , etc. Primes must also 

be placed on the parameters in the variational equations: 


1 d^’ _ B W 
k' dt Be' 


_ H. 

k' dt 


^ 1 .. 3W 

sin j' BN' 2 ^ Be' 


since they refer to the motion of the disturbed body. After differentiation the 
primes may be dropped without fear of confusion. In fact, primes are not needed 
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on X, j, since these quantities are not differentiated in the above equations; 

nor are primes needed on k and ; hence primes on them may be dropped be- 
fore differentiation. 

Terms in W which remain periodic after differentiation and after the primes 
have been dropped may be deleted at once from W, since we are interested in 
only secular changes in the orbital parameters. 

To illustrate the procedure the term 


2X'Y'3t ?H- 2 



- ^2 
2 


appearing in W is (Darwin 1880, Equation 37) 

[2(d'- 5)- 2fi] + 4F 772 ^2zr'2K'2 

t2(ff-- 5)+ 2£j]. + ^ |fi77^7t''‘ e-'^ [2(5--^)- 2fi] 

+ e>^2f + F^k^k'* [2(0*- 5) + 2fj]| 


Periodic terms have been deleted. If X and x' had been included in the 
above expression 2 (x - x') would appear in the e 2 q>onentials of the first term in 
curly brackets and -2 (x - x') in the exponentials of the second curly bracketed 
term; but since primes are unnecessary on x, these terms disappear. Also, 

& = 0 t + e and 8' = D't + e but D = D', so that - 8 becomes e' - e. 

We will now apply one of the variational equations to the terms in W in 
which the lag angle 2f ^ appears: 
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1 d<^ 

To find this term's contribution to 7- we must find 


B W 


2f, 


. Now e' 


k dt Be’ 

does not appear in tt'\ tz:'*, or ZI' so that the derivative operates only on 

the exponential terms . We obtain after differentiation 


3 W, 


2fi 


2 


Fi [2(€--o-2fi] _ ^ 4^.4 [2 (£'-€)- 2 fj] 


Dropping the primes, we have 


3 W 


2 ft 


3e' 


~ Fj sin 2 f j 7T^ 


where 

7T = Pp - Qq E ~ Pp " Qq 

and 


g2/M ft - e-2>^^i 

sin 2f, = ; 

^ 2 >/^ 


N is equal to zero when the earth and moon are the only two bodies in 
existence [as discussed earlier] , so that 77 = tt = cos (i + j)j . Also, in the 
case of a viscous earth Fj = cos 2ft, making the contribution of the Wjf ^ term 


. 1 d^ 

to 3T 
k dt 


1 

2 



sin 


4 




The other terms in W may be evaluated in similar fashion to finally give 
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= 

k dt 


(II-5) 


1 r 

^ ~ I w® sin 4 f j sin 4fj + Att^ sin 2 

- 4t 7^ sin 2 gj - 6tt^ K* sin 4h 
1 

j) 


as the secular rate of change of i, where k = sin (i + j 
Similarly 


dj _ h -At 

T dl - — [2 « sin 4 f, 


+ 77® /<r® sin 4 f + tt k"^ sin 4 fj 


+ 77 ® K® ( 77 ^ - K^) sin 4 h - -I- 77 ® /f [ 77 ^ - 3 sin 2 gj^ 


"2 77 K ( 77 ^ ” sin 2g+-277 K®(377^- k^) sin 2 gjj 


(II-6) 


gives the secular rate of change of j. 

Equations (II-5) and II-6) are respectively Equations (73) and (71) of Darwin 
(1880). 

The secular rates of change of the earth's angular velocity n and the inclina- 
tion i of the earth's equator to the invariable plane can be derived from (II-5) 
and (n-6) by application of the law of conservation of angular momentum: 

77 ® sin 4 f J + 2tt* sin 4 f + -^ k® sin 4 fj + 77 ® sin 2 gj 

+ 77 ^ ( 77 ^ - sin 2 g + 77^ K® sin 2 gj 

(n-7) 
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77 sin 4 fj 

+ -^ 7 T® /< (tt^ + 3 K^) sin 2 gj - 77 /< (t 7 ^ - /<^)^ sin 2 g 

1 3 

- 77 /<® (377^ + sin 2 gj - ^ sin 4h (H-8) 

These two equations are the last two equations of § 11 of Darwin (1880). 
Equations (II-5) through (II-8) are central to our discussion. 


4 77^ /<• sin 4 f ^ - TT^ (jt"^ - /<2) sin 4 f 
fl , 



CHAPTER m 


TIDAL FRICTION AND LARGE VISCOSITY 

A. Inclination of the Lunar Orbit to the Earth’s Equator 

In this chapter we discuss the secular motion of the moon for small 
inclinations of the lunar orbit to the earth's equatorial plane for earth-moon 
distances less than 10 earth radii. Equations (II-5) - (II-8) are not valid beyond 
10 earth radii because solar influence would have to be considered (see Chap- 
ter I, Section C). The effects of the sun beyond 10 earth radii are considered 
later in Chapter IV . 

The viscosity of the earth will be assumed constant throughout this discus- 
sion. Variable viscosity is considered in Section D of this chapter. 

Let us first write Equations (II-5) - (II-8) in slightly altered form. 

From Chapter II the moon's orbital angular momentum is 

^ 

The earth's rotational angular momentum is 

Lj. = C n 

The constant k introduced on page 21 is 

C _ 

^ GMm ^ ^0 

We wish to express our results with reference to a particular earth-moon 
distance c^; hence we can rewrite the first equation as 


27 
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-i 


M + m 


Mm 



bf 


where 




MX Mm Cj^ • 
M T m 0 


Since 


k 


— — Q r 
G M m 0 ^0 



But by Kepler's third law 


02 c 3 = G(M + m) 


Hence 


bk = C 


Equations (II-5) - (II-8) may now be written as 

* 

d Ljj ^ ^2 

~ ~2 ^ 4 f j sin 4f2 + 4 t7®/c2 2 

- 4 77^ K® sin 2 gj - 6 77^ sin 4 h] 


an-1) 
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dt 



TT^ K sin 4 f j sin 4 f + ^ 77^ /<^ sin 4 f ^ 


+ 



/<^) sin 4 h - K [tt^ - 3 t<^] sin 2 gj 


+ 


1 

■2 TT K (JT 


2 


sin 2g+ 47 tk®(3tt^-/<^) sin 2 gj 


(in-2) 


d Lg 
dt 



77® sin 4 f + 


2 77^ sin 4 f + 


1. 

2 


K® sin 4 fj 


(III-3) 


+ 77® sin 2 gj +77^ (77^ - sin 2 g + k® sin 2 gj 


di 

dt 


i! £ 

a Lg 


1 

2 


77^ K sin 4 fj 


77^ /c® (77^ - K®) sin 4 f - ^ 77 AT^ sin 4 f j 


+ 



K (77^ + 3 k^) sin 2 gj 



77 7 c (77^ - sin 2 g 


1 

2 77 T< 


5 


(3t7^ 


+ sin 2 gj - 



sin 4 h 


<in-4) 


Let ^ be the inclination of the lunar orbit to the earth's equatorial plane; 


then 


ip - i j 


77 - 


K - 



0 

cos ~2 


. 0 

sin ~2 
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For small (on the order of a few degrees or less) 77 ^ 1 and k « 1, 
Assuming small 4> and keeping terms through Equations (IH-1) - (HI-4) 
become 


dL„ \ ^ 1 

-JT - T “g~ JT2 ^ 


(m-5) 


dj ^ "^0 ^ 1 1 fl 7 . ^ 1 , . ^ ^ 1 c . ^ ' 

t = C — [2 ^ ^ 77 ^ /< sin 4 f j - 2 77 sin 2 gj (III-6) 

dLg j Pj 1 

— C 77 ® sin 4 f j + 77 ® >c^ sin 2 + 77 * sin 2 gJ (III-7) 

di ^ ^ 1 1 [1 7 . ^ 1 7 . r 1 7 . O 1 /TXT 

= ~Y~ — [^2 77 ^ X sin 2 gj + ^ 77 ^ X sin 4 f j ^ 77 ^ x sin 2 gJ (ni-8) 


where we have explicitly written for t. 


Usii^ the approximations 


- . 'P 4> 

X - sin ~2 ^ ~2 


'P ^ 1 'P^ 

77 - cos ^ % 1 g" 


We write Equations (III-5) - (in-8) as 


dt " 2 g ^1 


“ [( 1 “ sin 4 f j + i/;2 sin 2 gjj (III-9) 


2 

^ ^ ^ C ~ [sin 2 gj - sin 4 f 1 - sin 2 g] (ffl-lO) 
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dLg Tjj 1 pi ^p2 Jj2 

IT = C — [2 (I - sin 4 f, + ^ sin 2 sin 2 gj (m-11) 

2 

- T ^ ^ ^ 2gi + sin 4f, - sin 2 g] 4> (IH-12) 

neglecting powers of <p higher than 2. 

If Equations (III-IO) and (in-12) are added together we obtain the rate of 
change of 4> time: 

di/> ^ ^ ^ 

dt " 4 0 ^12 

(m-13) 




sin 2 gj 


Note that 


d^p 

dt “ 


'P 


d0 

so that for i/r = 0, ^ = 0. 

If the moon orbits the earth exactly in the equatorial plane, then the inclina- 
tion will remain zero. 

If the moon is slightly perturbed out of the equatorial plane so that i/r > 0, 
then the moon will move toward or away from the equatorial plane depending 
upon whether 


It * t) " t) 
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is a negative or positive quantity, i.e. 
d'/> 

(a) ^ <0 if (ni-14) is negative 

(b) ^ >0 if (III- 14) is positive 

For case (a) an equatorial orbit would be stable since small perturbations 
in would drive the moon back toward the equatorial plane. In case (b) an 
equatorial orbit would be unstable, because small perturbations in 'p would 
cause the inclination to grow at a rate proportional to i/». It is this second case 
in which we are mainly interested; we therefore want to examine (III-14) in de- 
tail to learn whether the tides can drive the moon away from the equatorial 
plane. 

We start with the coefficients of the sines of the lag angles, 
is always positive 

is positive for c < 21 earth radii 

Thus both these terms are positive in the region of interest. We next turn 
our attention to the lag angle terms. 

Equations (111-9) - (III-13) indicate that the tides which govern the evolution 
of the earth-moon system for small inclinations are the tides with speeds n - 2fl, 
2 (n - fi), and n, with the lag angles being gj, 2f, and g, respectively. These 
tides are called O, Mj, and Kj in Darwin (1883). 

To learn something of the nature of these tides we refer to Figures 6 and 7 . 
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Figure 6 shows ^ and n as a function of earth-moon distance for '/' = 0. 
Here use is made of the equations 

yc (M + m) 

q3/2 

which is Kepler's third law, and 

Lt " ^ 2 Lg L^, cos 0 

which is derived from the conservation of angular momentum. From this latter 
equation we obtain (remembering 0 = 0) 

_ Le Lj ■ Lm Mm 

“ "c" " c ^ c 

Figure 7 shows the angular speeds of the tides as a function of distance. 
The region to the left of the dashed line is inside the Roche limit (2.89 earth 
radii) where the moon would be torn apart by the tidal stress if the moon lacked 
cohesiveness; thus distances greatly inside the Roche limit are not physically 
realistic. Note that both n and 2 (n - fi) are positive for distances greater than 
the Roche limit, while n - 2Q changes sign at 3.83 earth radii (dotted line). 

The distance where n - 20 = 0 makes a convenient reference distance (at 
this distance the earth's rotation period is about 5.25 hours and the moon's 
orbital period about 10.5 hours). We henceforward take c^j as the earth-moon 
distance where n = 20: 


3.83 • a 
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where a is the present radius of the earth (6.37 x 10® cm). All quantities with 
zeros as subscripts refer to their values at this distance. 

We now examine the sines of the lag angles. Using the identity 


1 + tan^ 6 


and the assumption of viscosity we have 


sin 4 f j 


4 (n - n> C 

1 + 4 (n - n)2 


(in-15) 


sin 2 g 


2n C 

1 + nn^ 


(ni-16) 


2 (n - 2 0) C 

1 + (n - 2 0)2 ^2 


(m-17) 


where 


19 V 

2 g a/7 


and use has been made of the tangent formula for the lag angles (see page 20). 

The signs of (ni-15) - (III-17) are of the same signs as the respective 
speeds; thus sin 4f^ and sin 2g are positive while sin 2gj^ is negative for n < 20 
and positive for n > 2 0 . 

From the above considerations we can assert that 

d0 

< 0 for c < Co 


for all values of viscosity since each of the three terms in (111-14) is negative. 
Hence an equatorial orbit is stable at least up to 3.83 earth radii distance. 
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The sign of (III-14) for e > c,j depends upon the viscosity of the earth. To 
prove this, we examine this expression in the limit of low viscosity and then in 
the limit of high viscosity. 

In the limit of low viscosity 


where from Figure 7 


(n - n) 4 « 1 


(n “ il) S; lO'"* sec. 

This implies v « 10^® poises. 

In this case Equations (in-15) - (III-17) can be written as 


sin 4 f , 


4 (n - n) ^ 

1 + 4 (n - 11)2 ^2 


4 (n - H) ^ 


sin 2 g 


2n i 


1 + n2 ^2 


2n C 


sin 2gj = 


2 (n - 2 fl) C 
I + (n - 212)2 ^2 


2 (n - 2 n) c 


It is convenient at this point to introduce Darwin's notation (Darwin 1880) 


n 


\ decreases monotonically as the earth-moon distance increases (still remem- 
bering = 0) and obviously has value 0.5 at c^^; see Figure 8. 

Using this notation we have 
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sin 4fj %n^4(l-\) 
sin 2 g % n ^ - 2 


sin 2gj n ^ 2(1 - 2k) 


and (ni-14) becomes 


n I 


2 



( 1 - 


2X)- 4 



(1 - - 2 



This expression (see Figure 9(a)) was found to be negative by numerically 
computing the expression in square brackets for various distances. Thus, for 
small viscosities < 0 everywhere outside the Roche limit and an equatorial 
orbit is stable. 

In the limit of large viscosity for which 

(n - n> ^ » 1 


which implies 


» 10^5 poi 


poises 


Equations (III-15) - (III-17) become 


sin 4 f 


4(n - a) i 


^ 1 + 4 (n - Q)2 ^ V ~ ^ 
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sin 2 g 


2n i 

1 + 



2 (n - 2 fi) C 2 _ 1 ^ 

2gi = (n- 2fi) C ■ nC 1- 2\ 


This last expression holds only in the regions away from c^; see Figure 9(b). 
earlier sin 2gj = 0 at Cq = c while the expression above approaches infinity as c 
approaches Cq. The behavior of sin 2 gj near Cq will be examined later. 

For large viscosity (III-14) becomes 


1 


l± 


1 

Cn 

W 1 


Le/ 

1 - \ 




2 


The expression in square brackets is positive for c > c^; see Figure 9(b). 

We next inquire about the behavior of sin 2g^ near Cjj where n = 2 H. We 
make use of 


Let 


n 


]/G (M + m) 

c3/2 



C - C - Cq + Co - Co + X 


where 


X - c - Co . 

X measures deviations in distance from Co- The expressions for and n 


become 
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n = 


yC (M + m) 


fin 


, 3/2 


/ \ 3/2 . . 3/2 

(■•a ■•t 


Lt ^ C| 


•• (' • s 


- b 


(>• 4 ) 


n 


If — « 1, then by Taylor series 


3 

fi = fi„ - f ^ X 


Lj b 2 0^^^ 


n = 


Hn - 


0 2 Co C 


keeping only first order terms in x. 
So we have 


3fi„ 


n~2fi _ Hq 2coC *~2f^o + 


2 fio + 


b 

L ^0 ^ 


where 


0 + 


X 

n 


e = 


i 



b 


2 Cn 
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Equation (in-17) becomes 



Note that this expression is antisymmetric in x. 

If € « Cq (v » 10^5 poises), then sin 2gj has the features shown in Fig- 
ure 10. Sin 2 gi ranges from -1 at x = - e to 0 at x = 0 to +1 at x = +e . The 
peaks become sharper as the viscosity increases {and e decreases). 

Both sin 4f^ and sin 2g are only slowly varying for v » 10^® poises and 
are virtually constant between x = -e and x = + e. Further, both sin 4f j « 1 
and sin 2g « 1 for v » 10^® poises. Expression (III-14) is then seen to be 
zero for c slightly greater than c^, and the zero approaches c^j as the viscosity 
increases; we may then speak of (III-14) as being zero for c = c^ and positive 
for c > C(j with negligible error for large viscosities (»10^® poises). 

We conclude that becomes strongly negative for c < c^ and strongly 
positive for c > Cq. Thus for viscosities » 10^® poises an equatorial orbit is 
unstable for c > c^ and the moon will be driven away from the equatorial plane 
if perturbed. 

The transition of (111-14) from negative to positive for c > c^ was found to 
occur at a viscosity between 10 and 10^® poises by numerical computation. 

We summarize the major results of this section. 

(i) £ 0 for c < c^j for all viscosities and an equatorial orbit is stable. 

(ii) ^ £ 0 for c > c^ for viscosities less than about 10^® poises and an 

equatorial orbit is stable. 
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(iii) ^ - 0 for c > c„ for viscosities greater than 10^® poises and an 
equatorial orbit is unstable. 

These results hold for small values of yfi (the inclination of the lunar orbit 
to the earth's equator). 


B. Variation in the Earth-Moon Distance 

Let us turn our attention to Equation (10-9) and write it as 

1 '^ 0 ^ 1 
2 TV ^ P ■ 

(m- 19 ) 

© 

to show explicitly that we are discussing essentially the variation of the earth- 
moon distance in time and not the orbital angular momentum. 

In the limit of extremely small angles the terms can be neglected and 
(III-19) becomes 


- 

dt ~ 
^ = 


^ = 1 ^ 1 

dt 2 g b ^12 


(m-20) 


so that only one tide governs the variation in distance. If v « 10^® poises, then 

by the approximations of the previous section ^ « v; and if v » 10^® poises, 

then ^ . 

at 

The right side of Equation (III -20) is positive because sin 41^ is positive; 
therefore the moon is driven away from the earth. Also, is greatest when 
sin 4fj = 1, which occurs at a viscosity of about 10^® poises. 

Equation (III-20) can be integrated to give 
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1 Tq ft 2 


I3 - in = 2 75 


sin 4 f j dt 


If sin 4f j is only slowly varying, we have approximately 


1 1 "^o 

_L 13 _ ± 

13 ^^2 > - 7 a 


2 gb si" (tj - tj) 


for the dependence of distance on time. 

Returnii^ to Equation (HI -20), if the equation is divided into Equation (III-IO) 
and into (III-12) we obtain equations which eliminate the time: 


11 = i ^ _ sin 2g1 

d^ 2 Ly sin 4 f j sin 4 


di _ 1 b [sin 2g, ^ sin2gl 

d^ 2 Lg [sin 4 fj i sin 4f J 


If each side of these two equations is multiplied by kn, then 


11 = 1 ^ ^ _ _ sin 2 g 

^ d^ 2 Ljj sin 4 f j sin 4 f j 


di _ 1 [sin 2gi 2g"| , 

" d^ 2 sin 4 f J ^ sin 4 f j 


where we have used kb = C and Lg = Cn. 

Now because \p = i + j is very small, examination of Figure 5 shows that 


L-e i ^ Ljj j 


giving 
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0 = i + j = 
Substituting, we have 


dj _ 

i 1 


sin 

2 gi 

- 1 - 

sin 

2g 

d^ 

2 ' 

V 

sin 

-- 


sin 


di 

1 j 

{ '-'V 
' ^1 

sin 

2 El 

+ 1 - 

sin 

2g 

W ■ 

2 ' 

sin 


sin 




or finally 


kn 


d log j 
d^ 


1 


sin 2gi 

1 m 

sin 2 g 

2 1 


sin 4 f 

sin 4 f J 


kn 


d log i 
d^ 



sin 2 gj 

1 + — ; a~F~ 

Sin 4 f J 


sin 2 g 
sin 4 f J 


Using our previous approximations for sin 2g^, sin 4f^, and sin 2g in the limit 
of low viscosity we obtain 



These equations are found in §19 of Darwin's 1880 paper (Darwin 1880 pg. 312). 


In the limit of large viscosity 
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These last equations are found in §20 of the 1880 paper (Darwin 1880, pg. 317). 

Equations (III -21) and (III-22) will be discussed when analyzing Darwin’s 
theory of the moon's origin. 

We now return to Equation (in-19) and write it as 


- 

dt 


2 a b 



- sin 4 f j sin 2 gj] 


It is not generally true that 

sin 2gj| « I sin 4fjj 

for small \p for viscosities greater than 10^® poises because I sin 2gil may be 
on the order of 1. 

An example will illustrate this. Take 

V - 10^® poises 

sin 2 gj = 1 

ip = 3° = 0.0052 radians 

n - D = 1.66 X 10-* sec-l 

Then 

sin4fj % (n - D) I ~ 
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i/»^ sin 2 gj = 0 . 0027 
'//2sin4fj = 0.000006 


Obviously in this case 


sin 2 gj > sin 4 f ^ 


The i/t 2 sin 4 f ^ term is generally quite small and may be neglected giving 


if - 

dt “■ 


_1 

2 




0b^ 


12 


[^sin 4 f j sin 2gjj 


The point of this example is that even for small <p (on the order of a degree) 
neglect of the sin 2 g j term may lead to serious error. In fact, this term 
may have profound effects on the Ixmar orbit. We demonstrate this by examining 
some possible histories of the lunar orbit. 

Figures 11, 12, and 13, show (/;; sin 4fj,</;2 sin 2g, and sin 2gj; and 
respectively as functions of x for a viscosity of 10 poises. (All computations 
for Figures 11-14 were carried out with the computer program described in the 
next section.) The initial conditions are chosen to be i// % 0.4° at x = -0.4 x 10-3 
earth radii; it is labelled A in Figure 13. 

Since 


sin 4 f j + \jj^ sin 2 gj > 0 

for the chosen starting condition, > 0 at A and the radius of the moon’s 
orbit increases, so that the moon moves away from the earth (to the right in the 
figures) . As the moon moves toward point B, its outward rate of motion becomes 
slower and slower as sin 2g| becomes more and more negative. Past point B 
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increases and reaches a maximum at c. From c onward the rate of motion 
dt 

decreases. Note that ^ decreases for x < 0 and increases for x > 0. 

Now if the radius of the moon's orbit is initially less than Cq - e and the 
radius of the orbit is expanding, it must be that 



0 


at all points along its outward journey if the moon is to reach the outer regions 
pastCQ. In other words, 

sin 4 f j sin 2 > 0 

Now at X = -€, sin 2g^ = -1 and the above condition becomes 

sin 4 fj - ^ 0 

or 

(// ^ -j/ sin 4 f j 


at X = - e. 

Hence for initial conditions for which the radius of the lunar orbit is less 

than Cq - € and expanding, the above restriction must hold: the inclination 

must decrease below a certain value at x = -e for the moon to gain the outer 

d ^ 

regions. Thus, if such an orbit has a large inclination, ^ becomes small as 
the moon approaches x = - e and the moon "waits" near x = - e until tp has 
decreased enough to allow the moon past x = -e is negative for x < 0) 
and into the outer regions. 
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The net effect is that the distance x = -e acts as a "barrier" and will not 
let the moon through until '/' has dropped below a critical angle, which we label 

'Pc-- 



Critical angles for various viscosities are given in Table 2. 

Information regarding the history of the inclination for orbits which have 
ip > for X < -e is lost at x = -e , since <p must be less than or equal to 
in all cases to get past the barrier. 

Figure 14 gives an example of xp initially so large that 



and the moon moves toward the earth. It continues to do so until 

sin 4 f j + \p^ sin 2 gj = 0 

at point D in the figure. Here ^ chaises sign and the moon moves away from 
the earth. Thereafter the moon’s possible motion is as described before. In 
this particular case \p is quite large as the moon approaches the barrier, so 
that the moon must wait until \p drops below 2.7“ (the critical angle for 10^® 
poises) before moving into the region past x = 0.0, 

The effect of the barrier and the moon’s orbit shrinking and then e}q)anding 
can occur only if the radius of the lunar orbit is less than c^. If the moon 
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formed or somehow arrived at a distance greater than , then there is no 
restriction on the inclination and the moon moves continually outwards, since 
sin 2 gi > 0 for c > Cq . 

We summarize the major results of this section. 


(i) 


(ii) 


« - 1 r ■ 4,1 


for equatorial orbits 


^ g b 


[sin 4 f j + sin 2 gj] 


for viscosities greater than 10^^ poises and \p on the order of a 
degree 


(iii) \p < ^sin 4 f^ at c = ^ 

for an expanding orbit 

(iv) if sin 4f^ + sin 2gj < 0 (for c < Cq) an orbit will contract, and 
then expand. 


C . Computational Results 

The integration of Equations (in-5) - (m-8) was carried out numerically 
with a FORTRAN computer program using double precision variables.* The pro- 
gram is given in Appendix C. It was run on an IBM 360/91 computer at the 
Goddard Space Flight Center, as well as on a Univac 1108 at the University of 
Maryland. 

*This program was also used to obtain results similar to those of Gerstenkorn (Alfven 1963). 
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The program has the capability of integrating the full Equations (in-l) - 
(in-4), but the contribution of the neglected terms was found to be insignificant 
in the computations discussed below. 

The viscosity v, time interval At, and the initial values of t, i//, and ^ are 
read into the program. From the initial data i, j, n, and are computed, as 
well as Ly and L^. 

The program then iterates equations (in-5) - (ni-8) by computing the 

changes in i, j, L , and L„ according to the simple formula 
£ 



where X is 1, j, L^, or Lj^. The new values become 

^New “ ^Old ^ X 

is then recomputed from the new values and the process is repeated. At 
each step the new values of n, 0, fl, i, and j are printed, as well as Ai, Aj, 
and A\/r. 

After a certain chosen number of steps NQ At is adjusted so that the step 
change in ^ A ^ is constant for the remainder of the run. The reason for switching 
from constant At to constant A ^ is to insure that the time intervals at the begin- 
ning of the run are small enough so that the peaks in ^ near ^ = 1.0 are not 
missed (most of the runs start near ^ = 1.0). Later, as ^ increases and the 
change in ^ and the change in angles i, j, and in time become small, constant A^ 
is used to keep the run from becoming extremely long. 

If at any step |A i/i j or t A j/j 1 exceeds some chosen fraction called CRIT 


in the program the time interval for that step is halved and the step is repeated 
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until both 1 A i/i | and I Aj/j | are less than CRIT. The purpose of introducing 
GRIT is to avoid large changes in angle at any one step which would lead to 
cumulative errors after many steps. 

When ^ exceeds some chosen value XIMAX, or the total number of steps 
exceeds NLAST, the run is terminated. At the end of each run the total angular 
momentum is computed from the values of the last step of the run and is com- 
pared to the initial angular momentum. This serves as a check on how well the 
approximations used in writing (in-5) - (III-8) conserve angular momentum. 

After a run is completed its accuracy can be checked by halving the time 
interval of each step, doubling the number of steps, and repeating the run. 

The program also has the capability of integrating backward in time as well 
as forward. 

The program was run for various viscosities for which the moon is per- 
turbed from an equatorial orbit near c = c<j . The relevant data for these runs 
is summarized in Table 3. All runs stopped when the moon reached 10 earth 
radii distance from the earth; beyond 10 earth radii solar influence must be 
taken into account. No viscosities above 10^^ poises were considered because 
of the unrealistically long time scales involved. 

Figure 15 shows as a function of earth-moon distance for an initial per- 
turbation of 3° for viscosities of 10^® , 10^®, and 10^^ poises at c = c^. Note 
that \p decreases as a function of distance for 10^® poises, but increases for 
10^* and 10^^ poises. This behavior is to be expected from the discussion 
given in Section A of this chapter. 

The program was run next for \p-\p^a.tc = CQ-€ for 10^* , 10^^ , 10^® 
and 10^^ poises (Figure 16). This is the largest possible value can have near 



c = Co .without suffering some further perturbation from an equatorial orbit. 
Smaller initial angles at c = Co - e invariably gave smaller final angles at 
10 earth radii. 

Finally, Figure 17 shows i// as a function of distance for viscosities of 10^* 
and 10^^ poises for initial perturbations of 1“, 2“, and 3“. Figures 18 and 19 
give i and j respectively for the given initial perturbations. Curves for viscosi- 
ties between 10^® and 10^^ poises fall between the respective curves given in 
the figures. 

Note that only when the initial perturbation is about 3* does reach near 
10® at 10 earth radii as required in Goldreich's model; or equivalently, does j 
reach 6® at 10 earth radii. 

D. Variable Viscosity 

It was next assumed that the viscosity was not constant, but that the vis- 
cosity V was a function of absolute temperature T and that the earth was cooling 
down from an initially molten state. The purpose in doing this was to see 
whether the earth could cool off enough to be solid and have a high viscosity by 
the time the moon moved from the Roche limit (2.89 earth radii) to Cq (3.83 
earth radii). If so, the mechanism for driving the moon out of the earth’s 
equatorial plane may have been operative. 

The dependence at v on T was assumed to have the form 

V = (m-23) 

where 

Vq = a constant 

E* = activation energy per atom 
k = the Boltzmann constant. 
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A theoretical derivation of this equation is given by Glas stone, Laidler, and 
E 3 n:ing (1941) . We have ignored the dependence of Vg on T and have assumed it 
to be a constant. Experimental data shows that this equation holds fairly well 
for silicate melts (Clark 1966), with 

Vq % 10"'* poises 
E* ~ 2-5eV/ ato m 


Data on molten rocks are uncertain; the activation energy E* has approximately 
the range given above, but Vg may vary by orders of magnitude. 

A cooling law for the earth was required to give the temperature T as a 
function of time t. The law adopted here is derived in Appendix B. From Equa- 
tion (B-5) of Appendix B we take the form of the cooling law as 


T(t) 


Tg 

[1+ ZS^g3(t - tg)]"''' 


where 

Tg = temperature at time tg 

2 _ surface temperature 

average temperature of the earth 

12 77 O' 

® M C 

p 

where 

a = radius of the earth = 6.37 x 10® cm 

cr = Stefan-Boltzmann constant - 5.72 x 10“® erg cm“^ sec~^ deg"'* 
M = mass of the earth = 5.99 x 10^’ g 
Cp was taken to be 1.0 x 10^ erg/g-deg, giving 



52 


S = 1.46 X 10-20 

sec 

The computer program given in Appendix D was used to give the moon’s 
distance, and the earth's temperature and viscosity, as a function of time. The 
program differs from the program of Appendix C only in that the viscosity is 
allowed to vary in time rather than remain constant. The moon was initially at 
the Roche limit and in the equatorial plane of the earth, with the earth at a tem- 
perature Tq, Various values of E*, Z, and T^ were used to see if the earth 
could cool down near the melting point of rocks (about 1500“K) by the time the 
moon reached c q. 

The results may be briefly summarized. For E* 5 eV and > 10 
poises the temperature of the earth at Cq did not drop below about 1500" for 
initial temperatures between 2000°K and 3000°K with Z ^ to ^ . For E* % 

O A 

4.3 eV and Vq in the neighborhood of 10“^ poises, the temperature of the earth 
could fall below 1250" for the same ranges of initial temperature and values of Z, 
Apparently large values of E* and , which increase the viscosity at any given 
temperature , hasten the moon past c^ before the temperature has a chance to 
fall very low. 

Due to the wide variation in results and ignorance of the interval condition 
of the earth, it appears that we can make no definite statement as to whether the 
earth could cool down sufficiently from a molten state to have a large viscosity 


when the moon reaches c^. 



CHAPTER IV 


SOLAR INFLUENCE 


The history of the lunar orbit for distances greater than 10 earth radii, 
where solar influence must be considered, will now be investigated. Our discus- 
sion will be restricted to the behavior of J, the angle between the plane of the 
lunar orbit and Its proper plane. The angle j, the angle between the plane of the 
lunar orbit and the invariable plane, is essentially J for distances less than 
10 earth radii (see Chapter I, Section C). Hence at 10 earth radii we will join 
our previous solutions for j as a function of distance to those we obtain for J 
as a function of distance. 

Darwin obtains the rate of change of J with respect to ^ , with solar tides 
included, in Section HI of his 1880 paper. It was found by assuming the inclina- 
tions of the earth's equator and the lunar orbit to the ecliptic are small and 
applying the variations of parameters technique in solving differential equations. 
After a quite lengthy analysis he obtains (Darwin 1880, eq. 250, pg. 297): 


dlogj _ J_ 
kn 


if + a I 

•<*(/<, + a) (a' - /3' ) - a' b 7 — - b a | 

(K, I ' J 


(IV-1) 


]. (r (/<2 + a) + A + a) + bG - a D} 


where 


a - 


r' 1 

m + — 


T 2 \ e 


a = m /? = 1 + 


b = 1 
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Mq = mass of the sun 
Cq = earth-sun distance 

The change of i with time is given by (Darwin 1880, eq. 227, pg. 293); 


1 1 

k di ~ 2 g sin 4fi 


(IV-2) 


This is just the same equation as (III-20), where the only tide-raising body 
was the moon. The two are the same because solar tides and the direct gravita- 
tional force of the sun on the moon produce no secular change in the moon's 
distance. 

The computer program of Appendix E integrates Equation (FV-l) from 10 to 
60 earth radii for any desired viscosity. It assumes a constant step size in ^ . 

The angular velocity of the earth n is computed by assuming the total angu- 
lar momentum of the earth-moon system is conserved and that the moon re- 
mains in the equatorial plane of the earth. (The neglect of the frictional effects 
of solar tides and inclination leads to only small corrections in the final results.) 
These assumptions make the right side of the equation independent of angle, so 
that the solution of the equation has the form 


j, = . 1 ^ 

where F (v, ^) is the right side of the equation and Jj is J at ^ 

the initial distance A graph of J versus earth -moon distance c has the same 
shape for a given viscosity regardless of the initial value of J; larger or small 
initial values of J merely shift the curve up or down. 

The program was first run for the present-day values J = 5"9' and ^ = 3.96 
(60 earth radii) for a viscosity of 10^2 poises to obtain the small viscosity limit. 
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The result is shown in Figure 20 (dotted line). Note the close agreement with 
Goldreich's curve (dashed line), where the lag angles are assumed to be equal. 
The program was run again for a viscosity of 10^° poises. It showed negligible 
difference in results from 10^^ poises, so that the dashed line indeed repre- 
sents the small viscosity limit. 

The program was rim next for a viscosity of 10^® poises to extend the curve 
for j shown in Figure 19 for an initial perturbation in 4 > of 3®. The resulting 
curve is the upper solid line shown in Figure 20. The program was run again 
for a viscosity of 10^^ poises; it produced little change in the shape of the curve 
from 10 to 60 earth radii; hence the curve shows the large viscosity limit in that 
region. Note that the character of the large viscosity curve is quite different 
from that of the small viscosity curve. 

If the eaxth behaved as though it had a large viscosity from the time the 
moon was at 3.83 earth radii to the present time, then an initial perturbation in 
of about 2.5° at 3.83 earth radii would be required to give the present value of J 
of 5°9*. This is shown as the lower solid curve in Figure 20. However, viscosi- 
ties greater than about 10^’ poises give time scales of the orbital evolution 
greater than the age of the solar system. 

What is more likely is that the earth behaved like a liquid with high viscosity 
in its early history and then like an anelastic solid or liquid with low viscosity 
later on, which is what is observed today; so that the inclination J in Figure 20 
started on the upper solid curve at 3.83 earth radii and switched over to the 
dotted line, possibly somewhere in the region where the two curves merge beyond 
15 earth radii. Darwin (1880, §32, pg. 363) discusses the possibility of this kind 


of behavior. 



CHAPTER V 


DISCUSSION 

A, Critique of Assumptions Made 

We shall now examine the important assumptions made in obtaining our 
results for strong tidal friction. 

One important assumption we have made is that the orbit of the moon re- 
mains circular throughout its history, i.e. the eccentricity e of the lunar orbit is 
zero. The work of Darwin <1880, Section V), Singer (1968), and MacDonald (1964) 
shows that weak tidal friction decreases the eccentricity as we look back into the 
past until the moon reaches about 3 earth radii from the earth, where the eccen- 
tricity undergoes rapid changes. Since the present value of the eccentricity is 
0.055, this would imply that neglect of the eccentricity when the moon was at the 
reference distance of 3.83 earth radii would lead to negligible error in con- 
siderii^ weak tidal friction. 

However, use of Darwin's treatment of the eccentricity for viscosities 
greater than 10^’ poises indicates that e increases with time imtil the moon 
reaches about 16 earth radii distance from the earth; at larger distances the 
eccentricity rapidly decreases. This indicates that the eccentricity could have 
been large for earth -moon distances of less than 16 earth radii. However, a 
nearly circular orbit for the moon over its whole history is by no means ex- 
cluded. The earth could have behaved as though it had a large viscosity when the 
moon was less than 16 earth radii from it; beyond 16 earth radii the earth could 
have behaved as though it had a small viscosity. If this were the case, then if the 
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moon were in a nearly circular orbit at 3.83 earth radii, the eccentricity would 
slowly grow to its present value as the moon moved outward to its present 
distance. 

Another assumption which we have made is that harmonics higher than the 
second degree may be neglected in the tidal disturbing function (Equation II-4) . 

To show that this approximation is a good one, we note that the second degree 
harmonics in the disturbing function are multiplied by » where a is the 
radius of the earth and r is the earth-moon distance; this may be seen from 
Equations (H-3) and (II-4). If third degree harmonics were included in the dis- 
turbing function, then they would be multiplied by ; likewise, fourth degree 

/a\io 

harmonics would be multiplied by (— j , etc. Hence the third degree terms are 

reduced by a factor of from the second degree terms. For r = 3.83a, 

2 

= 0.Q68. Also, the contribution of the third degree tides to the rate of 
change in time of the inclination is small compared to that of the second degree 
tides (see the discussion in the last paragraph of this section). Thus the restric- 
tion to the second degree terms in the disturbing fimction leads to only a small 
error. 

We have further assumed that the moment of inertia of the earth C had its 
present-day value of 8.11 x lO"*® g-Cm^ = 0.33 Ma^. This implies that the 
earth's core had already formed. Darwin assumed that the earth was homoge- 
neous (C = 0.4 Ma^), as well as incompressible, etc., for reasons of tractability 
in solving for the response of the earth to the tidal force. Changing the moment 
of inertia to its value for a homogeneous earth would lead to only slight correc- 
tions in our results. 

Heating of the earth by the dissipation associated with the friction does not 
appear to be significant. The energy deposited in the earth as the moon moves 
from the Roche limit at 2.89 earth radii to 3.83 earth radii amounts to 1.43 x 
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10^® ergs. Assuming a specific heat of 10^ ergs/g-deg, the average change in 
temperature of each gram of matter in the earth is only 24“K. 

Our most crucial assumption was that the earth behaved like a highly 
viscous liquid {v » 10^® poises). Whether the earth could have behaved in this 
manner when the moon was at 3.83 earth radii depended upon the rheological 
properties of the earth at that time; these properties are unknown. 

O'Keefe (1972) points out that since the tidal potential varies like the in- 
verse cube of distance (Equation A-6, Appendix A), the tidal forces acting on the 
earth were 4000 times greater when the moon was at 3.83 earth radii than they 
are today, so that the material in the earth may have been near the elastic limit. 
In such circumstances the eartii may have behaved like a highly viscous liquid. 

At the present time the mantle of the earth responds to the tidal forces like 
an anelastic solid, with the tidal lag angles being small (MacDonald 1964). How- 
ever, the mantle responds to deformations of the earth's surface caused by ice 
loads as though it had a viscosity of about 10^^ poises (Gutenberg 1959, Chap- 
ter 9), requiring thousands of years to rebound after the removal of the loads. 
(This may be explainable in terms of diffusion creep; see Kaula 1968, pgs. 101- 
104). Now the period of the O tide with speed n - 2fi is given by 2 Tr/{n - 2D), so 
that the period ranges from infinity to about 5 hours as the moon moves through 
3.83 earth radii distance. Hence if the earth has a characteristic response time 
between these two extremes, then it should be excited by the O tidal force as the 
moon passes throi^h 3.83 earth radii. If the dissipation were great, then the IdLg 
angle of the tide would be large. Hence it is by no means clear that the earth 
would not respond as we have assumed, even with the present internal conditions 
in the earth, where the characteristic response time is thousands of years. 
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Our last important assumption was that tbe moon may have been perturbed 
out of an equatorial orbit by 2.5 to 3° at 3.83 earth radii distance from the earth, 
thus explaining the present inclination of the lunar orbit to the ecliptic. Whether 
the moon could suffer such a perturbation is not clear; conditions at that time 
may have been chaotic enough to produce it. However, several sources of the 
perturbation may be ruled out. The first obvious source of perturbation is a 
collision of a large meteoritic object with the moon. If such a collision occurred, 
then large amounts of meteoritic nickel might be expected to spatter over the 
moon’s surface*. Large amounts of nickel are not observed in lunar samples. 
The third degree harmonic in the earth’s figure will not give rise to long period 
perturbations in the inclination if the moon's orbit is circular. Further, it may 
be shown from Equation 38 of Kaula (1964) that the tides associated with the third 
degree harmonics in the tidal disturbing function give « 0, just as in the 
second degree harmonics, but that these terms are much less important than 
those discussed here. Also, the disturbance in the inclination caused by the pre- 
cession of the earth’s axis and the moon’s orbit may be shown to be quite small 
(«1°). The question of the source of the perturbation remains open. 

B. Relation of the Results to Theories of the Moon's Origin 

We will now examine how our results relate to the theories of the origin of 
the moon. The three principal theories, namely fission, accretion, and capture 
are reviewed by Kaula (1971). 

Darwin (1880) proposed that a primitive body rotating with a period close to 
its natural oscillation period was disrupted into the earth and moon by resonance 
oscillations induced in it by the sun. (That this was not at all likely was shown 


*I am indebted to Dr. O'Keefe on this point. 
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by Jeffreys 1930.) The moon would necessarily be thrown out in the equatorial 
plane of the earth. Darwin was forced to assume that the primitive earth had a 
very high viscosity to solve the inclination problem. He derived Equations (in-21) 
and {in-22) which give rate of change of inclination with distance in the limit of 
infinite viscosity. After commenting on the absurdity implied by the equations 
that the rate of change of angle was infinite when the earth rotated twice as fast 
as the moon revolved, he assumed that the viscosity merely had to be very large 
to increase the inclination of the moon's orbit to the equatorial plane from an 
infinitesimal disturbance to an appreciable angle. Darwin took this as the solu- 
tion to the inclination problem and let the matter rest. 

Our detailed analysis (Chapter III) shows that the initial perturbation in the 
inclination of the moon's orbit to the earth's equator must be about 2.5-3. O'* to 
explain the present inclination, with the viscosity of the earth being greater than 
10^^ poises. 

O'Keefe (1969) in his version of the fission theory suggested that the primi- 
tive body had greater mass and twice the angular momentum of the present 
earth-moon system. As the primitive body spun up, its figure progressed along 
the sequence of the well-known Jacobi ellipsoids and pear-shaped figures (Jeans 
1961) until it fissioned into the earth and moon. The system then lost mass and 
angular momentum through intense heating. While taking over Darwin's results, 
O'Keefe further suggested that even if the earth were molten after the moon and 
earth separated, the moon's orbital evolution would not begin until the earth 
cooled off appreciably, so that the moon would not arrive at 3.83 earth radii 
until the earth's viscosity was quite high. 

In Chapter III, Section D, we investigated the orbital evolution of the moon 
as the earth cooled off for a number of different activation energies and coeffi- 
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dents in Equation (in-23). In view of the wide variety of results obtained in the 
temperature of the earth when the moon arrives at 3.83 earth radii, it appears 
that the self -regulating mechanism proposed by O'Keefe does not exist. 

The accretion theory states that the moon formed from a ring of particles 
in orbit about the earth. The particles collided with each other and stuck to- 
gether, ultimately building up into the moon. 

The ring of particles would be expected to lie in the proper plane. The 
orbits of particles Inclined to the proper plane would process, thus lowering the 
chances of collision; at least all the orbits would intersect the proper plane, 
favoring accretion there. If the moon accreted from the ring much beyond 3.83 
earth radii, then the moon would tend to remain near the proper plane, so that 
its present inclination to the ecliptic could not be explained. However, if the 
moon formed in the proper plane within 3,83 earth radii (essentially in the equa- 
torial plane) , then the mechanism proposed here for driving the moon out of the 
earth's equatorial plane could have come into play. 

At any rate, regardless of how the moon arrived at 3.83 earth radii in the 
equatorial plane of the earth, whether by fission or accretion, if the viscosity of 
the earth was greater than 10^^ poises, and the moon suffered a 2.5-3.0° per- 
turbation in inclination at 3.83 earth radii, then the present inclination to the 
ecliptic could be explained. 

The capture theory has a simple answer to the inclination problem; the 
moon was captured in a highly inclined orbit to begin with and tidal friction has 
acted to decrease the inclination to its present value (Gerstenkorn 1969, 
MacDonald 1964). Of course, these theories begin with the present inclination 
of the moon’s orbit to the ecliptic and solve the equations of tidal friction back- 
ward in time to discover the/moon's inclination at capture. 
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These theories once again assume weak tidal friction with small lag angles. 
However, if the viscosity of the earth is greater than 10^® poises, and the moon 
arrives at a distance less than 3.83 earth radii in an inclined orbit, then the 
inclination must drop below the critical ai^le before the orbit can expand 
past 3.83 earth radii (Chapter III, Section B). Thus for large viscosities 
(greater than 10^^ poises) the orbit becomes nearly equatorial and we are faced 
with the same problem as before. 

C. Summary of the Important Results 

Assuming that the earth behaves like a viscous liquid in responding to the 
tidal force, and that the moon is in a circular orbit about the earth, with the in- 
clination of the lunar orbit to the earth's equator -S 20°, then: * 

(a) If the moon is less than 3.83 earth radii distance from the earth and 
the inclination of the orbit to the earth's equator is steep, then the orbit 
may contract and then expand, provided the viscosity of the earth is 
greater than poises. The orbit will expand monotonically if the 
viscosity is less than 10^* poises regardless of the inclination. 

(b) The inclination of the lunar orbit to the earth's equator will decrease 
or remain zero if the moon is closer to the earth than 3.83 earth radii, 
regardless of the viscosity. 

(c) The inclination of the lunar orbit to the earth's equator must be less 
than 2.7° when the moon is at 3.83 earth radii if the earth's viscosity is 
10^® poises; at higher viscosities the inclination must be even lower. 

(d) The lunar orbit will expand monotonically if the moon is at a distance 
greater than 3.83 earth radii from the earth, regardless of the viscosity. 

*<20^ $0 that Equations (lll>9) through (IH-12) hold. 



The inclination of the lunar orbit to the earth's equator will increase, 
or decrease (or remain zero) for earth-moon distances greater than 
3.83 earth radii depending upon whether the viscosity of the earth is 
greater than, or less than 10^® poises. 

If the viscosity of the earth is greater than or equal to 10^® poises, and 
the plane of the lunar orbit is perturbed about 2.5 to 3 degrees out of 
the equatorial plane of the earth when the moon is Just beyond 3.83 earth 
radii, then the present five degree inclination of the moon to the ecliptic 
may be explained. 



APPENDIX A 


DERIVATION OF THE TIDE-RAISING POTENTIAL 
AND TIDAL DISTURBING FUNCTION 

A. Derivation of the Tide-Raising Potential 

Consider the top diagram in Figure 21. The center of mass of the earth is 
located at point O; the center of mass of the moon is at Q; and the center of mass 
of the earth-moon system is located at point P. The earth and moon have 
masses M and m respectively. The vector h is directed from P to O. 

The center of mass of the system is taken to be at rest in inertial space, 
with the earth and moon orbiting about P in circular orbits. The angular velocity 
of either the earth or moon about P is CL. 

The X* y* z* coordinate system has its origin at O and is rigidly attached to 
the earth. The earth rotates about the Z* axis with angular velocity n with re- 
spect to inertial space. The vector r* = (x*, y*, z*) is the position vector of 
some unit mass element in the earth in the starred system. 

The bottom diagram in Figure 21 shows that the moon is located at r = 

(x, y, z) in the starred system, and the angle between r* and r is 0. s = r - r* 
is the vector directed from the mass element to the moon. We take |s | = s, 

I r* I = r* and | f | = r, so that the earth and moon are separated by a distance r. 

We wish to know the forces acting on the unit mass located at position r=^. 

— * 

Let us denote the total force on the unit mass as f.^. Let us further write f^ as 
the sum of two forces, one being the gravitational pull of the moon f , and the 
other being the sum of all other forces f (such as the earth's gravity, viscous 
forces, etc.). Hence we have 
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r. = f + c 

T HI 

represents the total force on the unit mass as seen from an inertial 

frame. 

We now wish to find the forces acting on the unit mass in the frame in which 
we have reduced the earth to rest, i.e. the forces as seen in the starred frame. 
Since the starred frame is non-inertial, fictitious forces will be introduced. 

Following Symon, Mechanics , Chapter 7, we write 

-*■ dn d^h 

^T* " f.j. - n X (n X r*) - 2n X V* - x r* — (A-2) 

Ij.* is the total force acting on the unit mass as seen from the starred 
system. The second term on the right-hand side of Equation (A-2) is the cen- 
trifugal force caused by the rotation of the earth on its axis. The third term is 
the Coriolis force, with v* being the velocity of the unit mass in the starred 
system. The fourth term arises from any variations in n. The last term arises 
from the earth's motion about the center of mass of the system. 

The third and fourth terms on the right-hand side of (A-2) will be assumed 
to be negligible, as they would be if the earth were changing its rotation rate 
only slowly and velocities relative to the earth were small. Equation (A-2) be- 
comes in that case 


ft* “ f + f-nx(nx r*) 


- 


d^h 

dt2 


where we have explicitly written f + f„, for fj* 


(A-3) 


Now 
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dh 

dt 


X h 


and 


d^h 

dt2 


fi x (O x h) . 


But n X (0 X h) = ^ where h = |h( = r. By Kepler's third 

law Q ^ gQ jjjay write 

r3 


d^h _ G(M + m) m r _ Gm - 
^^2 |.3 M + m r J.3 


d2 h 

We could have written this directly by recognizing that is just the 

acceleration of the earth's center of mass due to the gravitational pull of the 
moon. 

Equation (A-3) becomes 


r^* - r + - n X (n X >) - r (A-4) 

r j 

Let us now examine the f_^ term in Equation (A-4). L is the gravitational 
force of the moon on the unit mass located at r*. We can write f^ as the gradient 
of a potential: 


f_ _ V* V (x*,y*,z*, X, y, z). 


(A-5) 


Here 


V* = 




B X* 




3 ^ 3 


k* 


denotes the gradient operator operating in the starred system; i*, j*, and k*, are 
unit vectors along the x*, y*, and z* axes, respectively. V is the gravitational 
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potential of the moon at r*. Note that V is a function of both the coordinates of 
the unit mass and the coordinates of the moon, but that V* acts only on the 
starred coordinates. 

Taking the moon to be a point mass gives 

G m 

V = — 
s 

We now proceed to expand V in spherical harmonics about the center of the 
earth in the usual manner (Kaula 1968, Eq. 2.1.22): 

Since 

s = r - r*, s = (s • s)^ = (r • r + r* • r* - 2 r • r*)^ 

^ 'A 

~ (r^ + r*^ - 2 r r* cos 0)^ = r /l + - 2 cos 0^ 


We can thus write 



is of the form (1 + q)" where 

q " (^) - 2 (-^) cos 0 

and 

1 

n = - J- 
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Expanding (1 + q) " as a binominal series gives 


- „ . n (n - 1) ^ 

(1 + q)" = 1 + nq + q2 + 


so that our expression for V becomes 




2 r 


i ^ 1 

cos ® ~ 


where we have been careful to gather together terms in powers of — . This is 
nothing more than the familiar expression 


^ ■ “F" iL (cos ©) I— ) [ 

[m=0 ' ' J 


where (cos 0) is the Legendre polynomial of order m. 


If — « 1 we can ignore powers of — higher than 2 so that we can 


write 


V = ^ |i + cos 0 + I (-f) [cos20 - j]| 


Let us set 


G m 


_3 /r* 

2 V 




COS^ ® ~ y 


(A-6) 


which we will call the tide-raising potential. Note that is a second degree 


spherical solid harmonic function and 
harmonic. 


cos^ 0 - 


is a second degree surface 


We finally have 
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V = 


Gm Gm / r*\ 


cos 0 + V* 


If we again write f„ as the gradient of V we have 


f = V* 

m 


(^) COS0) 


+ V. 


The first term on the right-hand side of the above equation is zero because 
r^ = + y^ + z^ and nowhere contains x'*', y*, or z*. The function that the 

operator V* acts on in the second term can be written 


r* r* r cos © x*x + y*y f z*z 

cos 0 = = 


so that 


V* 


( x*x + y*y + 


z* z 


We are left with 


= Gm— + V* V, 

tn q t 


(A-7) 


Substituting the above equation in (A-4) gives 


f - f + G m — + V* - n X (n X r*) - G m 




or finally 


f^* - f + V* - n X (n X r*) 


(A-8) 
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Note that the first term on the right side of Equation (A -7) cancels the term 
associated with the motion of the earth about the center of mass of the system in 
{A -4), leaving the V*Vt term as the only term in (A -8) generated by the moon's 
gravity. We were therefore justified in calling the tide-raising potential. 

We note in passing that the centrifugal force may also be written as the 
gradient of a potential; 

- n x (n X r*) V* V^, (x*, y*, z*) 


where 

^ (x*2 + 

with n = Ini. 

In this case Equation (A-8) becomes 

= r+v*Vt + v*v^ (A-9) 

B. Derivation of the Disturbing Function 

In this section we show how deviations from sphericity of the earth give 
rise to a disturbing function. 

Let us again take the starred coordinate system to be fixed in the earth with 
its origin at the center of mass of the earth, and let S = (x, y, z) be the position 
vector of some mass element in the earth and A = (x', y' , z‘) be the position 
vector of some point E exterior to the earth (see Figure 22). Let S = |s | = 

(x^ + y^ + and A = |a| = (x'^ + y'2 + f = A - S is the vector 

directed from the mass element to the exterior point E and has length r = |ri . 
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The gravitational potential at (x , y , z ) is 


U (x', y', 2 * ) 




dV 


(A-10) 


where p is the density and dV the volume of the mass element, with the volume 
integral evaluated over the volume of the earth. 

The force on a unit mass at (x’, y', z' ) would be given by V U, where 


V = — — i* + i* + k* , 

Bx' By' ^ Bz' 


Note that the gradient operator acts only on the primed coordinates. 

We can e:!q)and r in a binomial series just as we did previously for s; 
Equation (A-10) then becomes 


U(x',y', z') 


= * (t) § (I) ■ t] * 


where '9 is the angle between S and A . 


If we neglect powers of I-t) higher than 2 we have 


TG p fG p ~ 

U(x'.y'.z') = J-A J A W «>s'PdV 


(A-11) 




dV 


The first term in Equation (A-11) is easily evaluated: 




GM 

A 
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This term gives the inverse square force. The second term of (A-11) is zero 
by virtue of having taken the origin of the coordinate system at the center of 
mass of the earth. The third term is the disturbing function R j; 


Ri (x', y', z' ) 


3 G r , 

1 

“ -ly — Ip 

COS^ - -r- 

2 a3 J 

L 3J 


dV 


(A-12) 


with the subscript I reminding us that V'Rj gives the force per unit mass in 
inertial space. 


Let us now write Equation (A-12) in spherical polar coordinates, with a and 
/3 being the longitude and colatitude respectively of the mass element, and a' 
and /S' being the corresponding longitude and colatitude of the exterior point; 
then 



cos^ ^ 

f r 

o 

o 

> 

o 



sin yS da dyS dS 


where 


cos ¥ - cos a' sin yS' cos a sin /3 + sin a' sin /3' sin a sin /? + cos /S' cos yS 


and 


dV - sin y6 da d/3 dS 


and Tj is the distance from the earth's center to the surface of the earth. 

Let us now write 

Tg (a, = a + a- (a, j3) 

where a is the mean radius of the earth and cr (a, /S) the "surface inequality" 
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which accounts for deviations of the earth's surface from sphericity, regardless 
of how those deviations arise. 

We introduce several assumptions at this point; first, that p is a function 
of radial distance S only and not a function of a and jB ; second, that a « a; 
and third, that a- {a, /3) may be written as a sum of second degree surface 
spherical harmonics fi) {t = 2). Of course any surface displacement in 

general may be expressed as a sum of surface spherical harmonics of all de- 
grees. We retain only the second degree harmonics since these are the most 
important. 

The second assumption allows us to write Rj (A, a , y3' ) as the sum of two 
terms , with the first term containing the volume integral evaluated over the 
spherical earth and the second term taking care of the "surface inequality": 

Ri(A,a’,;6') = p 8^ |cos2 - jJ sin is 82 da d]s d8 

^ Jo Jo Jq ^ 


.1 

2 


A3 


n 2'TT 


cr (a, yS) p 


COS 


2 ip _ 


3 


sin a2 da d^ 


where p^ is the density at the earth's surface. 

The first term vanishes by the first assumption because the integral over 
the angular part is zero. We are left with 


Ri (A,a',/3‘) 


3 

2 ^3 



COS^ V V 


r\j 

sin /3 da d/3 


Now by the third assumption cr (a, is a sum of second degree surface 


harmonics; cos^ ^ is also a sum of second degree harmonics. We then 
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recognize that the integral in the above equation is a sum of inner products of 
spherical harmonics. By use of the orthogonality of the Y„ (a , p )'s and keeping 
track of normalization constants we may finally write 

Ri(A,a',/3') = T^GPsa(f) (A-13) 

(See Kaula 1968, pgs. 65-67, for a general expansion in surface harmonics.) 

The disturbing function Rj at longitude a and colatitude /3‘ is seen to be 
proportional to the displacement of the surface at that same longitude and lati- 
tude; hence a body in the vicinity of the earth is acted upon by a disturbing func- 
tion which is proportional to the height of the displacement of the earth’s surface 
where the position vector from the center of the earth to the body pierces the 
earth's surface. 

If harmonics of degree n where n > 2 had been included in our expression 
for the surface displacement, then they would appear in our expression for the 
disturbing function correspondingly multiplied by For distances far 

from the earth « 1, so that these higher order terms are less important 
than the second degree terms. 

C. Conversion of Rj to R 

Equation (A-13) gives the disturbing function as seen inertial space. Gen- 
erally we want the disturbing function acting on the moon referred to the (ac- 
celerated) earth. We show how to find it below. 

Let Tj be the position vector of the earth in some inertial coordinate sys- 
tem. Likewise let r j be the position vector of the moon in this same coordinate 
system. (We still assumejthat (he earth and moon are the only two bodies in 
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existence.) Let = 1^2 " position of the moon as seen from the 

earth. 

If V is the potential of the earth, then by Newton's second law 


Mi?i = -mVV 

mrj - mW 


and 



M + m 
M~ 


^4 

vv 



is the acceleration of the moon as seen from the earth, and V 

represents the potential of the earth as seen from the moon. It is then clear 
that we must write 


R(A,a',y?') 


M + m 
M 


Ki 



(A-14) 


as the disturbing function acting on the moon as seen from the earth. 

D. The Tidal Disturbing Function 

The forces acting in Equation (A-9) displace mass on the earth; the dis- 
placed mass acts gravitationally on the moon and affects its motion. 

Assume that the earth responds separately to the centrifugal and the tide- 
raising force so that we may write 

cr (a, ^ = <7^ (a, (a, ^ 
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where cr^ is the displacement of the earth caused by the centrifugal force and a^ 
is the displacement of the surface caused by the tide-raising force. Equa- 
tion (A-14) then becomes 

R(A,a‘,;S*) = R^(A,a',/3') +R,(A,a',/3') 


where again the subscripts "c” and "f mean "centrifugal” and "tidal" respec- 
tively. We will call 





(f) 


(A-15) 


the tidal disturbing function. 

At this point we must take great pains to make clear the distinction between 
the tide-raising and tidally disturbed body. The two are not necessarily the 
same and must in any case be kept mathematically distinct to avoid incorrect 
derivations. We explain this below. 

Suppose we wished to find the action of the tides raised on Mars by Phobos 
on Mars' other moon Deimos. In the above discussion Phobos (the tide-raiser) 
is at point r = (x, y, z), and Deimos (the tidally disturbed body) is at A = 

(x', y’, z'). The force per unit mass on Deimos caused by Phobos' tides is 



3Rt ^ ^Rt 

i* -f 

Bx' By’ 


^ ^Rt 
j* + — k* . 
o z 


Even though depends on both x, y, z and x' , y z' the gradient operator acts 
only on the x' , y', z* coordinates. So much should be clear. 

Now suppose we have the case we are interested in, namely the action of 
the lunar tides raised on the earth on the moon itself. Here the tide-raiser is 
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also the tidally disturbed body, and the positions (x* , y z' ) and {x, y, z) are the 
same. We cannot drop the primes appearing in and £^jply the gradient 

B - 3 - 3 - 

_ — j* + — — 

3 X 3y 3 z 

to find the force per unit mass of the moon however, since the gradient operator 
can act only on the disturbed body's coordinates to retain the proper meaning of 
force per unit mass; thus the distinction between the tide-raiser and tidally dis- 
turbed body must be kept, even though they may be one and the same object. 

Darwin keeps the distinction clear in his 1880 paper by introducing the in- 
teresting artifice of giving the earth two satellites; the tide-raiser he calls 
Diana and the tidally disturbed body is the moon. When considering the action 
of lunar tides on the moon, Diana and the moon are, of course, the same object. 



APPENDIX B 


COOLING OF A PLANET BY RADIATIVE LOSS 

Let us suppose that a planet in empty space is cooling off by radiating heat 
into space from its surface according to the Stefan-Boltzmann law. Solar heating 
is neglected, and the planet is not surrounded by an atmosphere. 

Let us assume that the temperature distribution inside the planet has the 
form 

T(r, t) = T^(t)F(r) (B-1) 

where 

t = time 

r = radial distance 

T (r, t) = temperature at distance r and time t 
Tg (t) = surface temperature at time t 
F (r) = some function of radial distance 
Note that the temperature distribution has spherical symmetry. 

The planet radiates like a black body so that the amount of energy dQ given 
off in a time dt is 

dQ = - 477 RVT^^(t) dt (B-2) 

Here 

R = radius of planet 

cr = Stefan-Boltzmann constant 
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As the planet cools off, each element of mass dm inside the earth gives up 
an amount of heat 

dQdm = CpdT(r,t)dm 

in time dt, where Cp is the specific heat at constant pressure. 

The total amount of heat given off by the planet in time dt is then 

dQ = r Cp dT(r, t) dm 

^Mass of 
the planet 

This must be equal to Equation (B-2), so 

JCpdT(r, t)dm = - 4 t 7- R^cr T/(t) dt (B-3) 

Now from Equation (B-1) 

dT(r, t) = dT(t)F(r) 


gives the change in temperature with time, so that Equation (B-3) becomes 

JCpdT,(t)F(r) dm = dT^{t) JCpF(r)dm 

= -47TRVT;‘(t) dt 


Let 

J CpF(r) dm - I. 

Note that if we assume p has spherical symmetry, we may write dm = 
p{T) 4 dr to show dm as an explicit function of r. 



81 


The above equation can be written 


dT (t) 
T^t) 


This may be integrated to give 


- i [T,-=(t) - T-3(t„)] 



4 77 R^cr 

i 


dt 


4 77 R^a 

I (t " 


which may be written 


T,(t) 


Ts(to) 


1 + 


1277 R^O- 


Ts^(to) (t “ to) 


1/3 


(B-4) 


This expression gives the surface temperature as a function of time. 
If F (r) is known, the temperature at any point r at time t is given by 


T(r, t) 


T,(to) • F(r) 




77RVT/(to) 


( 


t - to)] 


1/3 


(B-5) 


Now suppose that Cp is a constant so that we may write 

I ” j'CpF(r)dm = Cp Jf (r) dm 

CpM jT,(t)F(r)dm 
" T (t) M 
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CpM Jt(f, t)dm 
TvCO M 


C 

p 


M 


<T (t)> 

T,(t) 


where M is the mass of the planet and 

fT(r, t) 

i> = ^ 


dm 


<T(t): 


M 


is the average temperature inside the planet weighted by mass. 
Set 


Z = 


_ surface temperature 
<T(t)> average temperature 


Z varies from ~1 to probably for any plausible temperature distribution. 

D 


Also set 


S = 


1277 R 2 (T 
MC_ 


S has a characteristic value for each planet. 
Equation (B-4) becomes using this notation 


T,(t) = 


Ts(to) 


[l+ ZS V(to) (t - to)] 


1/3 


(B-6) 


To give an example, for the earth 


R = 6.37x10* cm 


M = 5.98 X 1027 g 
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Cp ~ 1 X 10^ erg/g-deg 


S = 1.47 X 10- 


For Z = — and (t^) = 3000®K falls to half its value in about 1700 years. 
To demonstrate that a temperature distribution of the form 


T(r, t) = T,(t)F(r) 


(B-1) 


is not unrealistic, we note that the adiabatic temperature gradient is given by 


a = coefficient of compressibility 

g = gravitational acceleration 

The above equation may be rewritten as 


If the right side depends solely on r, we have 


(tJ " ■ Lc* 


(B-7) 


T = = T^FCr) 
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where 

= surface temperature 

and 

F(r) 



dr 


If the planet cools off in such a manner as to maintain the adiabatic tem- 
perature gradient at all times, becomes a function of time and 

T = T,(t)F(r) 

This is exactly the form which we assumed the temperature distribution 


had. 



APPENDIX C 


COMPUTER PROGRAM FOR CONSTANT VISCOSITY 

The computer program given in this appendix is discussed in Chapter HI, 
Section C, and in the comments listed in the program itself. 
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- * 

cc***^*****^* ******* 

|^"T40€rt 

C***** THIS PROGRAM IMEGRATES OARlllNS (13801 EQUATIONS TO GIVE THE 

C OF THE EARTH. 

C — T>^AR * l-N-4-1 880-^ 

C CN THE SECULAR CHANGES IN THE ELEMENTS OF THE DRBlT OF A SATELLITE 

Cr RE VOL V ING 8E0UT -A T tOAlA^ «ST0RT€O^ RM^ANET^ . ■ .. — - . — 

C IN 

- C— SCIENT SF^I^ T>APER8 -si R G ECRGg H O » ARP BARM IN » VG L^8^ PP^ - 

C CAMBR tOGE UNIVERSITY PRESS* 1908. 

C AND THE TVO PRECEEDING EQUATION 75*1 

-e - — — THE PAPER- CAN- A LSO B E FC U N O IN 

C PHILOSOPHICAL TRANSACTIONS OF THE ROYAL SOCIETY* VOL* 17(* I880i 

^ PP~^3 S8+* ^ 

C THE PROGRAM ALLOVS TWO APPROACHES - TO USE ALL THE TERRS IN 

C ORDER IN M«S1N(U4J>/2I* THE LATTER ME CALL THE SECOND ORDER 


C , THE EQUATIONS CAN BE INTEGRATED FORMARO OR BACKWARD IN TIME* 

- C --DEPENO^I N G U PO N W H6 TH 6 R MT I PE» 4-L OR 

C THE PROGRAM STARTS BY INTEGRATING WITH CONSTANT TIME INTERVALS 

C CONSTANT LM INTERVALS (CHOSEN IN THE PROGRAM TO GIVE CXls TO ASOUT 

— e ^yQOI-^- ALL T H I -S- A SS U M ES THAT THE PEP- CENT CHAN C E IN I O R ^ J I S 

C LESS THAN CRIT. IF THE FRACTIONAL CHANGE IS GREATER THAN CRIT* THE 

— € INTEIW^L--4S-H AL AM iP LP»T^i^THE-FltAC4i«NA^ IS-LESS^ TH AN- CRI T * 

C CRIT WAS INTRODUCED TO PREVENT LARGE CHANGES IN ANGLE TO AVOIC 

C COhSTANT OLM IS TO KEEP THE ITERATIONS FROM TAKING FOREVER* SINCE 


O F - OAR W IN VNOTAT-IOAU 

C****«C-2eR0 ^sREFERENCE DISTANCE 


C... ••SMALL KaC*(CMEGA-ZERCI*(C«ZEROI/(eiG G)8(B1G M)*(SMALL Ml 

C««... GOTHIC small Gss< 2/5 ( SMALL G)/( SMALL Al 

-■ C^*->^*B1C .GaUMlVERSAL GRAVITATIONAL XC.NSTANT 

C***«*SM^LL A^RADIUS OF THE EARTH 
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C 10**-4 /SEC. 

^ • V X & I nc^ Mrt 131*;^ oc^ wccrc t nc ■ irwff < rr^ a ^^nnc rww ■ r^B — ^ — * 

C PtANE OF THE I4O0NS ORBIT« PSi *I ♦ J. 

#* __♦ TC T Ut? iMi*l P n P TU F f *P TM *^ F m I ATTt D f * I f?l * ftiF Ahin TMP 

C INVAR I ABLE PLANE* 

C**«»*J IG T t C ANCL E-^- t ^ EE N- TliE-POeff^ O RB I T At^-PLAN€ -AND T HE-I N VARIABLE -- 
C PLANE* 

C»»»**P S I t l i A N D a APE I N nADlA N S * A fl O XJ AP€ TM€ S A A € A NGL£^- 

C IN DEGREES. 

- € ... * * 0P6I t Cl t AN D O J AR£ -T>4fe- €HA NG C S- I N T H S BBS P EC TIV - E ANGL ES * 

C#«««*VIS IS THE EARTHS VISCOSITY IN UNITS OP 104A16 CGS* 

C*». .. L X— IS THS-ORaiTALr ANGULAR Ng^EAtTW-QF^ THE S YSTE M- I N -UlilTS—OF- — — 
C IOA440 CGS. 

THE EARTH IN UNITS OF - - 


C I0AA40 CGS. LC*C*N. 

- CL H -ANO^-CLe-ARE THE RES PEC T1V€ -CHANG ES I N -L4<-ANO^-LC^ 

C LT IS THE TOTAL ANGULAR MOMENTUM OF THE SYSTEH IN UNITS OF 

— e TOAAAG^ ess-r — — — 

IS THE »*QHENT QF INERTIA OF THE EARTH IN UNITS OF 104*44 CGS. 

— e C ^ S M ALL - K4 B . — — — — 

C««-..T IS THE TIME IN UNITS OF 1C449 SEC* DELT IS THE CHANGE IN T. 

C4 4444 C N 0 SECTIO N . — 

C 

--€- y .-.-w^ ^A i TS - iS N ALL K» » < TAU - ZE R O>44a^< G OT W lC S HA LL G) U N QARWINS 

C NOTATION! TIMES B CAS DEFINED HEf£ ) IN UNITS OF 10**31 CGS. 

C.a...B IS S 0:>TM e iG G4S M ALL AI/C B IG Mi S H ALL MH *< B IG M l *C S MALL — M4 *D S I N 
C UNITS OF 10**40 CGS. LMs8*XI. 

<-w^w^ ir AE™ES™T9 j^< 24 C SMA LL G>»CS NA LL A I44S NA LL W)H N UN IT S OR 1 0* * >^ ia CGS. — 
C.....A3 IS SORTCCBIG G>*CB1G H4SMALL MI/CSHALL A)**3^/DS**3 IN UNITS OF 

— C L0**-4 SEC.— GNEGA«A3/^<-Xl**^T^ 

C 


C.....NRUN IS TFE NUMBER OF DATA CARDS TQ BE READ. ALL INITIAL DATA FOR 

— C A SI N GLE RU N I S C0N TA1W60 CN A SIN GL K C AR P. 

CRIT IS THE MAXIMUM CHANGE PERMITTED IN THE ABSOLUTE VALUES OF 

C r «« . . M TI M E— H FOR A -B A T C H O F R U N S - INTE GR AT CO F O R WAR D - I N- TI M E . AN D 

C MTIME*-1 FOR INTEGRATION BACKWARD IN TIME. 

C.....ANGLE - INITIAL VALUE FOR PSI»I * 3 IN RADIANS. 

C . . . .^C-«i--IHl-TTALr- V A LUE OF X C» 

VlSF s VISCOSITY OF THE EARTH IN UNITS OF 10**16 CGS. 

oeuTl F- ^ S T^ F 4N-T^M€— IN- UN lTS O F -M>*»0^- SBC.- — 

C DELT IF SHOULD ALWAYS BE POSITIVE* 

€ .**. *XI MAX » T H E M AXI MitN V A L US C F X I T Q WH IC H T H E ■ PRO GRAM -IN T EG R A TES ^ 

C THE RUN STOPS WHEN XI = XI WAX* 

C...«*NP s THE NUMBER OF ITERATIONS DESIRED USING CONSTANT STEP SIZE IN 

.^C TI ME. 

C..«..NF^1« NL* 1 MEANS THE PROGRAM USES THE SECOND ORDER APPROXIMATION. 

-■C-**-.-**MF»3. NL^2 MEANS THE PROGRAM USES ALL THE TEBIiS IH-OARWINS 

C EOUATICNS. 

G.**.*NF^1» . NLFtg MEANS THE PPOCRAIC BUNS BOTH TMF SECOND ORDER 

C APPROXIMATION ANO DARWINS FULL EQUATIONS ON THE INITIAL DATA* 

C APPROXIMATION BY HALVING THE STEP SIZES ANO DOUBLING THE NUMBER OF 

— C ST EPS A N O R EP CATI N G T H E RUN.-- IF T HE CHEC K IS NOT O B S | R E O. NC1*G. > 

C.*...NC2 PERFORMS THE SAME FUNCTION AS NCI FOR THE FULL EQUATIONS. 
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isr« 0002 


IS f 0003 - 

ts^ 0004 

ISh 0006 
IS f-0007 

— 000 6 - 
I‘S^ 0009 
— 0040 
ISK 0011 
-I SI^ 0042 
tS^ 0013 

OOt-4 

ISK 0015 
tS#^0016 
IS^ 0017 
-4S4 OOIO^ 
ISK 0019 
IS4 0020 
ISI^ 0021 
—464 00a«^- 

^_4S4 0022- 

- - 46 K 0020 

- 164 0025 - 
ISt*. 0026 
IS4 00 26 
ISK 0030 

I S4 0034- - 

IS^ 0032 

134-0033 

ISf 0034 
134- 003S - 
ISK 0037 
— ^-364-0038 
IS^ 0039 
IS 4 0040 

-134 0041 
ISK 0042 

134-0043 

tSK 0044 


_X3K 004S 
ISK 0046 


r» -yNL^lST^ 13- T H g Or Pg^Ml TTEO ONg WNr THE 

PUK TERMINATES IF NLAST IS EXCEEDED* 


C 

OQUatE PRECISION DZERO • C • L T *A I . A2 * A3 » A4» B * ANGLE • XI « VIS. LM »CC • 

1 — L-E-^N'rGM£ €A rP^*X * TF-4 1 rTG-rSG 41G1 *SG l-.OLrE# Dt4« L70*T rOELT-rOS-r 

2 XIMAX 

DOUBLE PRECISION AY^OI ♦03.0ELT4 - - 

DOUBLE PRECISION TF .SF *TF2 .SF2. TG2. SG2.TH. SH 

DOUBLE PRECIS ION^ PSl-rOP 54- - - 

DOUBLE PRECISION XlF ♦ VISF . DELTIF.TSTART 

— DOUBLE PREC ISIW T4-.T2^T3.34^TS*T3.^T7* T3«T 0*T4O*^T 14^732^X43* T 14* 

1 T 15.Tl6.T17.Tia. Tl9fT20*T51.T22.T23.T24,T25 

30U8LE-PREC434CN -SS*SJ*3*4 — 

DOLBLE PRECISION L .R1 «R2.P 3.D1 1 .0 J J 

COL3LE- PRECIS ION- L4 - - 

DZER0a3-fi33e730SDC 

3s6*4-4C0 - -- - — — - 

LT*34.200 

OS^SORT4C2ERC4 - - — 

Al*l*3 IS 157D44C/DZEB0446 

A2=2.T56DO - - — - 

e=3.681701D04DS 

—33^ V2 *49 1E5D0 /l OZEPOODS 1 — - 

U1«0*0C0 

^ CP I T^T IME - 

1 FORMAT ( I5.FI0*2. J5I 

OD -400 NR^ l .NRUN - - - 

C44444THIS SECTION READS IN THE INITIAL DATA* ANGLE IS IN RADIANS. 

- - READ <5*2 > ANGLE *4IF* VI SF* DEL44F^X1MAX,T3TART-*NP*4F-*NL*NC4*NC2* 

1 NLAST 

2 —4DRMAT<DO*2*04l*2-*4D40-a-* 15*412*174 

C44444ENC SECTION* 

CO - 1 00 -4 AstNF*-KL - 

IF (NA *£0« 14 ND^NCl 

IF XN4 *€Q* a4NO*NC2 . - - - 

NCs 1 + NO 

DC- 300 41B« 4.NC - 

NCKECKss-1 

40 . XlasXlF . - 

NC*NP 

3F 3410 ♦EXl-. a I NO»24KIP . 

VISsVISF 

OeLT4^DELT4F - — — -■ 

TsTSTKPT 

34»VIS4A2 - • - ^ 

C44444THXS SECTION WRITES OUT THE INPUT DATA. 

WAITE C6*fil VIS - - 

6 Format < IMI.///. lOX. lOHVISCCSlTT^.OlO. 2. lOH 10O416CGS.///) 

WRITE <6 *483 3^ELX3F-*34MAX*jyC*Jtf:^tilL *NC 1 *NC2^NLAST, — 

18 FORMAT (5X.54DELT1*.D10.3.1X.9H10449 SEC .5 X *6HXT MAX^ «C1 0 • 3 . 5X* 

4 3HN a= * 1 5 * 5 X . 3HNF « * 1 1 * 5 X • 3 HKL * ♦ X 1 * 5X* AHN C l - * 11 *.13.*5lXj 

2 6HStt_AST»,I5.///> 

WRITE <6.521 CRIT*MTIME_ - 

52 FORMAT < 5X. 5HCRIT».F1 0.4*5X.6HMTXME= » 1 3. /// I 

CROW *« END SECT ICM* . - — 

C44444ThIS SECTION WRITES THE HEADING* 
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IS^ 004a 


C *v* * **l,L OU4NTITICS W tSrvCN ttt ^rt1C-”#RTO1W4, 

C NOTE THAT UNDER HEADINGS PSI* 1, AND J XPSI« XI 1« AND XJ ARE 

-e P P I N T C G t gi ViN€ AL L -ANGbES 4>^-DeGR g €S * AL S O , T+4E-F ITTST^O^- AN0~ 

C DJ LISTED {READING FROM LEFT TO RIGHTS IS FOR THE INFINITE 

€ VISCChS MTf 4HM^T^ THe-+WFHI^TE -VlGC^SI Ty G H A N GE— TN--ANGLE- TS 0E 

C CQ^PAPED TO THE ACTUAL CHANGE IN ANGLE (SECOND 01 » OJ HEA0lN6«> 

7 FORMAT C6Xf 4HTIME *9X« 2HXI 1HN*6X«5H PSI «3X » 5H0MEGA t 7X «2HDI t 1 OX 


ISX 0049 
1S> ^050 
IS^ 0051 
tSN -0052 
ISN 00E3 
-tSH 0054 
ISh 0055 
— IS-N 0050- 
IS^ 0057 
—iOfr-OOM 
ISN 0059 
— 15 N 0050- 
XSN OOfil 

W ^ A. ^ ^ ^ 

£ 9 T* WXrmC 

XSN 0053 
- fSN 0 050- 
XSN 0065 


^N- 00 6 6 


ISN 0067 
-^NeOOOO 


C4****END SECTION. 

C * ** *4 T H45-SGC^T-IDH COMRWTEG THE INI T1 AL VALUES OF LE -i 
NT*0 

-L-M«-BOX-|- ■■ - - - - — ^ 

CC^OCOSC ANGLE 1 

PSI«ANGL6 — • ■ 


Nf ON e OA^ -l^^ AND ^ 


XPSIas l80*04ANGLE/3.14t5S 

15 — D€tT«0€fc7 4 - 

17 CONTINUE 

L€*-tM4CC-4 050RT4ffcl40CCI40e-4--LT442 - LN44624 

xiat.M/'e 

N»UEyC : — 


0MEGAsA3/<X|443> 

-5S«05-»K{-ANGLE4- 

S3«SSyCS0RT(SS442 N ICC 

3»DAR5 INiS34 

1*ANGLE-J 

^ W m A A AAA 


4 tM/LE)0«2} 


XJal80.00J/3« 14159 

-- COAHNOOENe-SeCTteN^ 

IS USED ONLY HERE TO FILL IN ZEROS* 

ItR ITE — <6^3 vH* HP SX * O ME GA -r 4 H-»4^W 

CC00N4N4N4N44N4 40044 

. A * A A AAA A A A A AA AA A 

^■^p-^^p’^r^p^p'*' — *"“ - . ■«— - »■ -• 

5 AY«0*5004PSI 


*-l-4-0e TO 2 6 


I SN -007^ — 
ISN 0071 

ISN 0072 
-4SN“007^3— 
ISN 0074 


C44444THIS SECTION COMPUTES ALL THE TERMS IN DAR4INS EQUATIONS* 
K«DSI5<Ay^H — - 


P=»CCO$(AY) 

fOMPUTg THE TANCFNTfia filN Rfi HF ^ljc i aJS. a er*» 

^ » "V w m r‘*^iF f 1^ w r tw 1 rf • O f O f W B * 16^ 9 " ^ WmsC Al # ■ 

TF«a*04N4A4 

- SF*a*C4TP/X W04TP-4Oa4 ^ — 


TF2ai2«04(N40M£GA)4A4 

t^04TF 2/4^04 T«0 4^- 


I5N 0076 

- 1 < A 

»3rl9 W T T 

ISN 0076 


ISN 0060 

y <6 ^i*> ^ ■ 

M~9~r V V U 6 

ISN 0062 
IS N 00 63 
ISN 0084 
X B N 0066 


TG2^( N 42 *04CM£GA I 4A4 

. 1 i*> *1 o » 

V 61 • •4F^f46AiP"6 

TH«2.040MEGA4A4 

04T44/4 047144624 

TF|32*04CN^CMEGA)4A4 

S A 4 ^ ^QA TF in/l 4^04TF6442T - 

TG«N4A4 
■fiC»2 ■ C41 

TGX«(NH-2*040ME6A|4A4 


ISN 0066 
ISN 0067 
“' IS N 0060 
ISN 0069 


C CONFUTE THE TERMS* 


T2»2*04P4444K4444SF 
Ta»0 *6 D04 N 44 6 4 S F2 


T 4 >P44 64X4*4 24 SGI 
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l$^ 0101 
t Sf» 0t0 2 
ISh 0103 


ISK OIOS 


tSf 0107 

K 4# V9W 

IS^ 0100 


ISh 0112 

I SN »i»3 

XS^ 0114 

ISh Oils 

ISI^ 0116 


1$^ 0119 


1S> 0121 


XS4 0123 


IS4 0125 


ISh 0127 


ISh 0125 


T20«P4*34K4*34<P*42-K**2)*SF 


T2 2»0 « 5004P44 €4K4 < P44 2^ 3* 0*K4*2 > OSS I 

- T 2 3 »0^ 00 0 *P4K4(< P44a^K » 42 i 44 3l » S0 — — 

T24»0 *500 «P4K«« 54 <3«04P«424K442 >4S62 


C44444EN0 SECT1GN» 


; FULL EQUATIONS. 


OLW=0.5004A1*CT7-T64T9-T10^TI t |40eLT/C XI ♦♦ 12) 

& 3* >> Al»< ' Ti a» TI3 -l> TI 4» Tl B-^ TI 6* T17 4 T1 8)* DgLT/<L M» Xl»412 i 

01 3kA 1 4 ( T 19-T20-T2 l4T22-T23-T24-T2S)4D€LT/€ LE4XI4*! 2) 

60 TO 27 

26 K*CSIN<AT) 


C««. • .COMPUTE THE TANGENTS* SINES OF THE LAO ANGLES* 

y Fl » 2*0 »<N-OM 60A)» A4 ^ 

SF1*2 .04TFI/1 1.0 4 TFI442) 


S6«2.C4T0/( 1.0 4 T0442> 


S61=2.04T61/€ 1.0 4 TG1442) 


SHs2.04TH/( 1«0 4 TH442) 


SF^2.C*TF/I 1.0 ♦ TF442I 

^ ftici I'T c t i-c '*caa*c - 


T1s0.5004P*454SF1 



ISh 0134 


ISh 0136 


ISK 0136 


ISK 0140 


T7^P663*K4434SF 


T9-0.5004P4474K4S01 


T 1 1«0 . 5D04P44 54k4 SG 


T 13-0 .50O4P4474K4^l 



91 


1S^ 0142 
tStr 0^>3- 
IS^ 0144 
014 5- 
1S^ 0146 


Tt 4=P4454K«434rSF — ““ — “““ 
TlE=0 •500*P4474K*SG1 

3*5^1 

T1 7==0 •5D04P4474K4SG 
-7ta*t-*^Dil4P4454K44»34S€^ — 
Tl9=l •5004P4434K4434SH 


ISK 0147 
iS h 0^140 
IS^ 0149 
-iS^0i06 — 

-15^ 019+ - 


C«4444TH1S SECTION COMPOTES THE CHANGES IN LE* 1_M* J* AND I FOS THE 

^ ^r%w%^i^r\ A ♦ 44 

^ w It 0VIV V ■ ■ 

36 DLE*-A 14 <T14T24T3 140ELT/C X I** 12) 


-0t M»414»0^ SD04T T 44T5 MiDE+T / C XI * 4 1 21- 


fSh 0152 
ISK 0153 
+0K 01 5 4 
IS4 0155 


DJ=A140*500*< P*47*K*SGl-P4*74K*SFl-P4*5*K*S6>40ELT/fLM*Xl4*l2) 
C44444END SECTION. 

2^ ^ - ■ - 

C44444THIS SECTION INSURES THAT OI/I CP DJ/J NEVEP EXCEEDS CBIT. 
0^ )- 19^«»30- - 


IF (RI - 0.0) 47t46«46 


IS)> 

iSf 

0156 

0157 

47— 

46 

Pl« -PI 

CONTINUE 


ISK 

0150 

0159 


IF <ftj - 0.0) 49«40.4e 


I9K 

IS4 

-0+60“ 

0161 

49 ■ 
46 

CONTINUE 


*19 h 
1SN 

0 10'2 
0163 

34 

Tr 1 — CR'IT f 34«34»35 

IF <RJ - CRIT) 19.19.3S 


ISA 0160 

35 

OEL T* OCL. T XS 1 U 

IF <NA .EO. 1) GO TO 36 


ISK 

V X W f 

0169 

19 

t W A — 'vCv. ~"K~W Ow ■ -fTS- O-.'" 

CONTINUE 




o c 1 a 

C44444TH1S SECTION HALVES THE STEP 

SIZE VHEN NCI OR NC2 EQUALS i* NCHECK 

ISK 

0170 

- — ^ 

“ fM«ur ^ 1 A ri^ • I 1 t pl. — v 

IF 4NB .EO. 2) GO TO 41 
GO TO 42 - - 


— f-^rA- 

ISK 

01T3 

41 

NCKECKs-NCHECK 

-ttrt- • 1 /i#i 


ISK 

0176 

^ « •’.’w 

4 9 

GO TO 42 


ISK 0176 


IF (NA .EG. 1) GO TO 36 

t C t KIX CA O i CO 3tfl 


ISK 

0102 

42 COHTINLE 




C44444THIS SECTION INCREMENTS THE 

IMPORTANT QUANTITIES. 

4 9 T% 

ISK 

0104 
0 196 


1*1 ♦ Cl 
J* J 4 C J 


ISK 

0106 


Xt 1*1 80^041/3. 141 59 

■w ft* 1 a A. A* .AtRC 


ISK 

0188 
0 1 B9 


PS1»P51 ♦ DI ♦ 03 

VOATx f tV'l^ lAlSO _ 


ISK 

0190 

A « A I 


LE*LE4DLE 


1$K 

0192 


N*LEyc 


ISN 

0194 


XX*LM/B 
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ISh 0197 

T ^ A ^ 

■■ S 9 

ISK 0199 

«> oa« 

I9^ 0201 

#1 _ 
1 ^ ^ V 

IS» 0203 
fSK ^204 - 
IS^ 0205 
-ISK^ € W OO 
IS^ 0207 


1S^ 0206 
.4.^1^.0209^ 


ISf 0210 
4 s» o a i I 
ISK 0212 


XSN 021$ 


ISf 0217 


l$^ 0219 


I$^ 0224 


XSI^ 0226 


tS^ 0229 


, — Tm^ T -t>0 €tT — — ^ 

C*#* •♦NT=HUM8ER CF ITERATIONS DONE SO FAR IN A RUN, 



C499**ENC SECTION, 

COOOOOTf^TS S^GTT Q N CO MPU TES T H E C^ANO ES IN— T AN&-^ IN T H E LI M IT -OF- 

C INFINITE VISCOSITY fOARUlN 1380 PAGE 317,1 

r^^^<ONPUT*T4EN-4~S-NOT eEGUN-UNTTk^-X4«4^000I TO AVOID DIVIDING BY Q, 

IF (XI - l.OOOlDOl 3It32,32 

31 O4'1«O^O0O — 

CJJ«0«ODO 

— EO-TO-23. ^ , 

32 t>C90ME6Ay<LT^MI 

R^l »4, 041^14 .0-H. l^ C-UO-2^,94t4 , 

R2-0.5C04 f 1 .0/CL.T-LM) 4 1,0/LHI9<1.0 4 Rl) 

R3*^0 ♦S € » 0 * < 1 ,0 /CkT - L NV 4 I , 0/L N I »<4^»O^14 - 

01XSOLMAR29T 

C JJaOL I4*R 3 4J - 

33 COKTXNtE 

.C44A44ENO SEC7«H^ - — - . 

C«**«*TH1$ SECTION PRINTS OUT THE N£M VALUES FOR THE IMPORTANT 

Q q ua n titi es^- — 

VRITE < 0 , 2 } T,XI ,N,XPSX tOMEGA »DtI ,DJJ ,Xi t ,X J,DI «DJ«OPSl 

3 AT ^4SX,^04 9 »4,-iX,01 3,6, I X, 01 0, 3., tX,F7,3, 2C4X»Oi0,5T,.lX«O10,4^> 

1 2UX« F9,3),3( IX,D10,4I) 

- C*966*EN0 6ECT40N^ — 

C96664>ThtS SECTION DECIDES WHETHER CONSTANT DELT CR CONSTANT OLM SHOULD 

C STATEMENT 12 GIVES CONSTANT STEF SIZE IN TIME, 

Gy^^^ TAT EM SNTfi 1 3 AN O $ 0 G IVE CO N S T A NT STE P 4N~ V«U — 

IF CNT • NO 112,12,13 


IF CNA ,EQ, 2} GO TO 50 

€0__W>.._i4 

50 OELTaKCABS<0,00720a9009Xl6612/(AI*0,5D09( TT-T84T9-T1 0-Tl 11 1 1 

—51 CCN TI N U E- 

GO TO 14 

12 OELT^OELT-1 

14 CCNTINCE 

-C66 9 66 EM C ^CT40K, - 

C SHOULD integration BE FORVAFD OR BACKNARO IN TIME? 

T R - ( MTIM E .,€0^- .*14.-OeL.T*»OELT 

C,«,«,1S NLAST EXCEEDED? 

IE-4NTwNLASTl 4,4 »99 , 

XIMAX EXCEEDED? 

— 4 l ^ -AX . 

C44444THJS SECTION COMPUTES THE TOTAL ANGULAR MOMENTUM AT THE END OF THE 

- C -TUIAU-^TT— SERVE S^S-A^CHECK-CN HON MEU THE I TFII AT TON _3CMF Mg MORkS- 

99 CC«DCCSiPSX» 

4^T C *D S O R T g. E»42 4 LM 6 6a4a,4>6.LEAL M6CCT 

MRITE C6,U) LT.LTC 

14 FO R MA T /J^liOX,lftHTN|TTAI. _ AKG* MnM,»,D14 ,6, 1 flX, 1 6HFINAL AKG, 

1,014,6) 

--499- C ON T INUE 

STOP 

END- 



APPENDIX D 


COMPUTER PROGRAM FOR VARIABLE VISCOSITY 


The computer program given in this appendix is discussed in Chapter HI, 
Section D, and in the comments listed in the program itself. 
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. 

CC*4t****«*4t«4e*«**^*# 

^CM*JMiTHlS.^lS^.JlOEZ- ' 

THIS P«aSftAM INTEGRATES DARWINS Cl8flC> EQUATIONS TO' GIVE THE 

-C -EVOI.UTION _DF__THE HQQNS f CTRCAII l^gl.jQ-RBlT FQ ? _ A_ VA R 1 ABLE VISCOSITY 

C 

-C— 

c 
_c 
c 

c . 
c 

.x.,_ 

,c 
. c 
c 
c 
c 

c 

X . 

c 

c 

c 
c 
c 

__c 

c 

_c 

c 

c 


OF THE earth. 

QARJ*1M_U86QI IS 

ON THE SECULAR CHANGES IN THE ELEMENTS OF THE OR8IT OF A SATELLITE 

REVOL VI NG AS QU T A T IDALLY DISTOR TED PLANET 

IN 

SCI ENT I FI C P_APERS^_BY . SIR GEORGE. HOWARD PARWIN* VO L. Z* PP 208":3R2 

CAMSRIDGE UNIVERSITY PRESS* 1908. 

JrHE^QUATlQNS^_AgE_J:m. ■PGS-J^AJ(Ll24J^ANO__2»^ CFQ UATIONS 7t * 73# 

AND THE TWO PRECEEDING EQUATION 75.) 

thK PAPER _CAN_ ALSO 8E FOUND. IN 

PHILOSOPHICAL TRANSACTIONS OF ^HE ROYAL SOCIETY# VOL. 171. 1630# 

PP 71.3 - 691.* .... : 

THE PROGRAM ALLOWS TWO APPROACHES - TO USE ALL THE TERMS IN 

DARWINS EQUAT_IOtiSjL_ OR KEEP ONLY TERMS UP TO AND INC LUDING SECOND 

ORDER IN K=SIN( (I>^JI/21 « THE LATTER WE CALL THE SECOND ORDER 

^ AHPJlQXimT LOK# 

THE EQUATIONS CAN BE INTEGRATED FORWARD OR BACKWARD IN TIME, 
■DEPENDING UP ON WHETHE R M TJJj£.^>L_OR -1* 


THE PROGRAM STARTS BY INTEGRATING WITH CONSTANT TIME INTERVALS 
IT DOE S THIS FOR NQ STE PS. AFTER NO STEPS IT SWITCHES OVER TO 
CONSTANT LM INTERVALS (CHCSEN IN THE PROGRAM TO GIVE DXI» TO ABOUT 

0*00 1*1 ALL THIS ASSUMES T HAT T HE PER CEN T CHANG E IN I OR J IS 

LESS THAN CRIT. IF THE FRACTIONAL CHANGE IS GREATER THAN CRIT* THF 
_JNTE8 y AL. I S H ALVE D UNTIL THE FRACTIONAL CHANGE I S LESS THAN CRIT. 

CRIT WAS INTRODUCED TO PREVENT LARGE CHANGES IN ANGLE TO AVOID 
CUMULATIVE ERRO R. THE S TRATAGE M OP SWITC HING FROM CONSTAN T PELT TO 
CONSTANI DLM is TO KEEP THE ITERATIONS FROM TAKING FOREVER# SINCE 
AT LARGE XI 0X1 I S _ QUITE S MALL F OR C QNS TANT PELT. 


C 

C...-*SQME OF DARWINS NOTATION. 

r - - 

--.r-zifDn — nTQTAMrp 


C.. ...CMEGA-ZERO-QMEGA AT C-ZERO 


C 9 • 

^^.FMAtL KsCM( omfga— ZFRO>* ( C— zppo) y (niG GI*(BTG miai'^mall 

Ml 

c«« 

•*.TAU«Z£P0-(3/2>*C8IG GI*(SHALL M> /(C-ZERO >**3 

C.J 

... gothic RMALL G— < >*( small GJ-Z^tSMALL Al 


— C*. g. 

...BIG G=UNIVERSAL GRAVITATIONAL CONSTANT 


c. . 
C«* 

••• SMALL AsRADTUS OF THE EARTH 
...BIG M=MASS OF THE EARTH 


C 99 

«« .small MsMASS of thf moON 


. Ca* 

...SMALL G*GRAVITATT0NAL CONSTANT AT THE EARTHS SURFACE 


r - * 

---W t «5 THF nPN.S fTV HP TMF FAPTM 


C.....WE HAVE SUBSTITUTED BIG G FOP DARWINS MU ABOVE. 


C 

C**^tATHIS SECTION DEFINES THE MOST IMPORTANT QUANTITIES. 

C.« 
c .. 

...XI IS SORT! EARTH-MOON DI STANCE/RFFERENCE DISTANCE). 
...OZERO IS THE REFERENCE DISTANCE# HERE IN UNITS OF EARTH 

RADII# 

c 

Ctg 

AND WHERE N-2FOMEGA. 
...OS-SQRT(DZEROI * C“ZER0=SMALL A*0S. 


c*^ 

T< TMF OnTATTniJAt AMGtII AO VFt nrTTV RF THF FADTM TKl 

c_ 

<I.E. MULTIPLY THE VALUE GIVEN IN ThE PROGRAM BY 10A9-A 

TO GET THF 

c 

VALUE IN CGS UNITS.) 
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O R B IT AL ANGULO Vg L OCITY OF T Hg WO OM t W UN IT S OF 

C 10«^4 /SEC* 

,.r.^ J»&I ■ I.S. TMF AMfil F flPTMFFM THF FAPTHA PniiATnr > TAI Pi AMF AMD THF 

C PLANE OF THE NOONS ORBIT* PSl^l + J* 

THE AN6LE BET WEEN THE- FAftTHS FQItA^nPlAI. Pl.ANF- AND JW£ 

C INVARIABLE PLANE* 

J 13 THE A NG L£ ~a£t M £ E N THE MnnWS ORBIT AI PLANE ANO THE INVAP I ABl E 

C PLANE. 

C *****. PSf* I * A NO .i ARE -. I N RA D IANS. XPST.XH* AND XJ APF THF SAME ANGLES 
C IN DEGREES* 

C*. **.nPST-»-nT* AND nj are the C HA MGFS in TMF PFSPECTIVg ANGLES* 

C VIS IS THE EARTHS VISCOSITY IN UNITS OF 10AA16 CGS. 

-C^*..**^L!l..LS-_TME__QftBlTAL ANGULAR MOMENTUM OF THE SYSTEM IN UNITS OF 

C I0AA40 CGS* 

r*****tF. IS THE POTATinNAI ANfilH AR, MOMENTUN OF THF FARTH IN UNITS OF 

C 10AP40 CGS* LEsC4N. 

X^*-»*-* OLN -Ald^OLE ARE THE RgSPgCTtVF CHANCFS TN I M AKiH I F. 

C«***«LT IS THE TOTAL ANGULAR NONENTUN OF THE SYSTEM IN UNITS OF 

C 10444Q. CGS* 

C**«**C IS THE MOMENT OF INERTIA OF THE EARTH IN UNITS OF 10*444 CGS* 

-C CaSMALI — KAB* 

C««***T IS THE TINE IN UNITS OF t0**9 SEC* CELT IS THE CHANGE IN T. 

CAAA44END- SECTION*- 

C 

C*.**^*A1-IS (SMAU K>4lTAt»»7Fani44PyrfinTHtr- SMALL Gl flN QAWMINS 

C NOTATION) TIMES 0 (AS DEFINED HERE) IN UNITS OF 104431 CSS* 

IS S QBTUBT G G ASMAi t AlViRTfl IH-SliAl I MI)4f BTG MI4(aMAi l Ml4nS IN 
C UNITS OF 104440 CGS* LM*84Xl* 

-C*^ * **A P I S i a/ (2 AfSM A U - OAI S M A H - AlAC S M AL L M)IIN UNITS OF 1044-17 CGS* 

C A3 IS SQRTUeiG G)4(8IG M4SMALL Ml/ (SMALL A)443)/DS443 IN UNITS OF 

X. 10AA^_C_SBC^-QMESAai A 3/< X 1 4 , 4 31 ■ 

C 

rAAAAATHE TIME VARIATIOM fiF THF VlSCflSlTY IS GIVEN BY 

C VI$=VlSZ4EXP(Be/TEMP) 

VIS7grnFFFTCTFNT OF VtSm SlTY IN UNITS QF 1QA41& CGS* 

C«* *«*eB3:AC7IVATI0H TEMPERATURE IN DEGREES KELVIN* 

-C^***TEMP IS THE TEMPERATURE QF THE EARTH IN DEGREES KELVIN* 

C*.*..THE TIME VARIATION CF TEMPERATURE OF THE EARTH IS GIVEN BY 

-C TEMP^TFMP2/I Cl *04BETA4.1TEMP2 4 »31 f sLT-^ Y START)) p 4Q^333) 

C** ♦**TEMPZ=TEMPERATURE OF THE EARTH AT TIME TSTARTi GIVEN IN DEGREES K* 
-C* * * * ♦ BET A«3A.iAAEA O F EARTH)4( STEFAN^BQLTZ* CQNSTI/( BIG M4SPECIFIC HEAT) 
C TIMES 

_C SURFACE TEMP*/AVG* T EMP* QF EARTH* ROUGHLY FROM 1/2 TO MAYBE 1/S* 

C FOR SPF* f^AT-10447 ERG/GRAM4DEG BETA»1 «47410A4><^1 1 /DEG44341 0440SEC 

C**** *Tim(FB IS THE LOWEST TEMPERATURE PERM ITTEO* IF TEMP DROPS BELOB 
C TLOVER IN THE EQUATION FOR TEMP* THEN TLOMER IS USED IN THE 

X. VISCOSITY EDUATTON^ TLQMgR IN DEGREES KELVIN* ^ 

C 

„CA* 44.4 THIS. SECTION EXPLAINS THE INITIAL INPUT DATA* 

C*.**.NRUN IS THE NUMBER OF RUNS TO BE MADE* ALL INITIAL DATA FOR A RUN 

X IS READ FROM TMQ CONSECUTIVE CARDS* 

C*«*«*CRIT IS THE MAXIMUM CHANGE PERMITTED IN THE ABSOLUTE VALUES QF 

-C DI/I_ AHD_DJ/>I IN A SINGLE STEP * 

C FOR A BATCH OF RUNS INTEGRATED FORWARD IN TIME. AND 

_C MT1ME«-1 FaR INTEGRATION BACKWARD IN TIME* 

C«*«**ANGLE ^ INITIAL VALUE FOR PSI^I IN RADIANS* 

C«.***X IF a INITIAL VALUE OF XI* 

C*«.**VISF VISC3SITY OF THE EARTH IN UNITS OF 104416 CGS* 
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C VlSF IS NQT USED- I N T H IS P WOGWA »l » S O A NY VALUg HA Y Bg U S E D * 

C.*.*.DELTIF = STEP SIZE IN TIME IN UNITS OF 10#*9 SEC. 

C- ^ DELT-IE- SHOULD- -Ai^MA.YS-fl£_PnS IT I VE . ^ 

C.....X1MAX = THE MAXIMUM VALUE OF XI TO WHICH THE PROGRAM INTEGRATES. 

C THE RUN STOP^ WHEN XI- =„ XIMAIL^ 

C.....TSTART * INITIAL VALUE OF TIME IN UNITS OF 10**9 SEC. 

C^. ITERATIONS DE SIRED USING CONSTANT STEP SIZE IN 


C TIME. 

^-CL.-.-^lfel^J4L=l^EAN5 T HE PR OGRAM USES TH E SECOND ORDER APP R OX IHATIQN* 
C.....NF=2* NL=2 MEANS THE PROGRAM USES ALL THE TERMS IN DARWINS 
C EQUATIONS. ■ 


C.....NFS1* NL=2 MEANS THE PROGRAM RUNS BOTH THE SECOND ORDER 


C 

C. 

X- 

c 




ISN 000 3 
TSM OflQA 
ISN 0005 
ISN 0QQ6 


ON THE INITIAL DATA, 

iNCl^l GIVES A CHECK ON THE DATA RUN FOR THE SECOND ORDER 
APPROXIM A TION B Y HAL VING THE STEP SI ZES_ AND OQUBLI NG NUMBER OF 

STEPS AND REPEATING THE RUN. IF THE CHECK IS NOT DESIRED. NC1=0. 

C PERFO RMS T HE SA ME FUNCTION AS NCI FOR THE FULL EQ UATIONS. 

IS THE TOTAL NUMBER OF STEPS PERMITTED IN ANY ONE RUN. THE 

C*4(*A«END SECTION. 

> DOUBLE PRECISION Q ZERO .C^XT.A1^A2 . A3 . A4 .8 « ANGLE.XI . VI S .LM , CC t 

1 LE.N.OMEGA.P.K.TF I.SFl .TG.SG.TGl.SGl .DLE.DLM. LTC. T. DElT f DS . 

2 X IMAX 

DOUBLE PRECISION AY.Dl .D J.DELTl 

^UBLE PR gCISlQN-JFfSEjT-F2iSF2,T G 2t SG2fTHfSH 

DOUBLE PRECISION PSI.OPSI 

QQUBLE PRECISION X IF . VlSF .OELTl F .TSTART 


ISN 0007 


ISN OOOB 
_-LSmiQD9. 
ISN 0010 


DOUBLE PRECISION T I .T2.T3 .TA.TS.Te.TT.T8.T9.TI O.Tl I.T12. T 1 3.T14. 

I _Ti 5 tT le » T.17 f >Tj£P-a.T. 2| j T22 jT2 3 . T25 

DOUBLE PRECISION SS.SJ^U.I 

lyXJBLg PRECISION L.RI.R2 .R3.0H .DJJ 

DOUBLE PRECISION TEMP.TEMPZ.8ETA .TLOWER.BB.El . EE. VISZ 


ISN 

ISN 

0012 

0013 

C^d.llDO 
„ . LT«34.2D0 

ISN 

0014 

DS^DSQRT(DZERO) 

ISN 

0015 

Al« t « 31 31 S7D44C/0ZER0446 

ISN 

001 e 

A2«2.756O0 

ISN 0017 

a=3«68t701004DS 

ISN 

0018 

A3^l 2 .491 8500/1 DZERQ4DS ) 

ISN 

0019 

Ul^O.OOO 

ISN 

0020 

READ (5.11 NRUN.CRIT.MTIME 

ISN 

0021 

1 FORMAT (lS.Ft0.2,I9> 

ISN 

0022 

00 100 NR*l.NRUN 

C44444THIS SECTION READS IN THE INITIAL DATA. ANGLE IS IN RADIANS. 

ISN 

002 3 

READ (5.21 ANGLE •XIF.VI5F.D6LT1F.X1MAX.TSTART.NP.NF.nl. NCI .NC2. 
1 NLAST 

ISN 

0024 

2 FORMAT C09.2fDll.2.4010«2.I5.4I2.I7l 

ISN 

0025 

READ (5.521 BB. VISZ.TEMPZ.BETA.TLOWER 

ISN 

0026 

52 FORMAT (5D10.41 

C44444END SECTION. 

18N 

0027 

DO 100 NA^NF.KL 

ISN 

0028 

IF INA .EO. 11 NO*NCl 

ISN 

0030 

IF (NA .EO. 21 NO-NC2 

ISN 

0032 

NC» t ♦ NO 

XSN 

0033 

DO 100 NB*l.NC 

ISM 

0034 

HCHBCF=»^I ^ 


ISN 0035 


XI»XIF 


40 



ISK t>03e -• - -NOSW- 

ISN 0037 IF (NB •EQ. 2) NQ^2*NP 

MX&aUlSF- 

ISN 0040 OELTI»OELT|F 

ISN 0041 T=TSXA« - ; 

C4**%4TMIS section COMPOTES THE INITIAL VISCOSITY, 

0042^ - ^ - - TEMPaT£MP2/C f L^^EXAiLi T£JtflZ*A31JLLT^TfiTAPT 1 I ♦♦Q, 33 3 1 . 

0043 IF (TEMP - TLOWERl 50,50,51 

0044 SO TEMP*TL^«ER 

0045 51 E1=8B/TEMP 

0046 - ..EE;?nEXRtHU - 

0047 VIS3SVISZ4EE 

C44444EN0 SECTIDN* 

ISN 0046 A4=VIS*A2 

C4*44r.ilLXHlS^ -SECT-JQM_ MJ &t TES HUT _THF TMPtlT P ATA^ 

ISN 0049 tfPiTE (6,6> 

ISN 0050 6 FORMAT -LLHiX ■ 

ISN 0051 WRITE (6,53> BB,V1SZ,BETA 

LSN- 00S2 S3. FORMAT LS3_>3HBRg,m O, A , I X , 7HDFG RFFS, 5X ,5HVTS7^,01 0 > A> 1 10^104416 C 

lGS,5X,5HBETAA«D10,4tlX, I 0M/104A9 SEC,/) 

I SM 0 0 5 3 TFMPy.T^ 

ISN 0054 54 FORMAT (5 X, 6HTEMPZ^,0 I 0, 4 • 1 X, 7HDEGft£ES,SX ,7HTLaVER *«D 10 ,4, I X,7HDEG 

- 1REES,/Z.) 

ISN 0055 WRITE C6,18> DELTi F,XIMAX ,NQ,NF , NL,NC1 ,NC2«NLAST 

ISN 0056 - TB FORMAT. C5 X ,6 HDELTl= ^ m0,3 »J.X^ 9 hLL Q . 4 49 SEC, 5X,6HXIMAXa, D1 0,3, SX, 

I 3HNO=, I5,5X,3MNF* ,1 1 ,5X, 3HNLs;, H •5X«4HNCls* II , 5X,4HNC2=t 1 1 ,5X* 

.3.„.6HhU,AST«,I5,///> 

ISN 0057 WRITE (6,57) CflIT.MTIME 

ISN OOSB- - 57 FORMAT- t6Y>gMrpiT« > Fi o_fc ^ Ry tmfst 1 

C44444END SECTION, 

-- C.44 44*THlS_BFCTTON _WaiTES THE MgAOTNO^ 

C ALL QUANTITIES PRINTED OUT ARF IN UNITS GIVEN IN THE PROGRAM. 

C 4IOtE-:EHA3^UNQRP HFAOTMGS PSl , 1, AMn J KO SI , XII, ANn KJ ABF 

C PRINTED, THUS GIVING ALL ANGLES IN DEGREES, 

ISN 0059 , WRITE -I5,71 

ISN 0060 7 format (6X« 4HT I ME, 9X.2HXI .OX, IHN ,BX «5H PS I , 3X iSHOMEGA, 7X, 3HVIS, 9X 

C44444END SECTION. 


ISN 

ISN 

-ISH. 

ISN 

ISR 

ISN 


CBSBBSTiil S,. SECTJIQN _C Q MPUIES_JriiE_ INITIAL VALUES OF LE, N, OMEGA. I, AMD J 

ISN 0061 NT=0 

— XSN. 0062 LM:=B4XI 

ISN 006 3 CCaOCOSCANGLE > 

ISN 0064 J5SlgANGLE_ 

ISN 0065 XPSI=180.04ANGLE/3, 14159 

— I SN 0066 IS OEl TaPELTl 

ISN 0067 17 CONTINUE 

— TSN 0066 LE^LMACC^ ^ _Q5QRTC C LM*CC144g ^ LTAAg ^ LM442I 

ISN 0069 XI»UM/B 

- ISN -0070 N??L£yC 

ISN 0071 0MEGA=A3/(X1*431 


ISN 0073 

SJaSS/DSQRTC SS442 ♦ (CC ♦ LM/LE > **2 1 

— 

ISN 0074 

J»0ARS1NCSJ1 


ISN 0075 

l^ANGLE-J 


ISM 0076 _ . 

XllJLl_BlL.a«lZ3mJ.*lS9 


ISN 0077 

a 

XJ>180«Q4J/3, 14IS9 


ISN 0078 

WRITE (6,31 T, Xt,N,XPSI , OMEGA, VIS, TEMP, xn,XJ 
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■ ISN 

T «N 

0Q7.R 

nn«n 

5 A^Y=sO.SDOAeSI - 

C***A#TH1S SECTION COMPUTES THE VISCOSITY. 

TEMP =TFI«a7yf f 1 ^Q*ttFTA*fTgMP7AA31AfT^T ST ART 1 1 ♦ * ft . 33 3 1 


ISN 

OOSl 


IE I TEMP - TLOWFRI 56.55. S6 


-tSAL 

nos? 


TPMPaTl nv^P 


ISN 

O0S3 

S6 

El*BB/TEMP 


T^W 

QflRA 


FFsDFiCPiFl) 


ISN 

0C85 


VIS=VIST*EE 




. CA**** END SECTION. _ . 


ISN 

0066 


AA*V1STA2 



I SN a087 lF_.fl*A mEQm. .1 1 , CLQ TQ . 

C***^*tmS SECTION COMPUTES AtL THE TERMS IN DARWINS EQUATIONS. 


T<tlU ftftRq 

- Krf)Sl HI AY 1 


XSN 0000 

P^DCOSCAYI 




COMPUTE XHE-TANOENTS. SINES OF JTME^J_A£_ ANGLES. 



ISN 0001 

TF=2.0*NAAA 


ISN 0092 

SFa2.0*TFyll.0fTF**2J _ 


ISN 0093 

TF2®2.0F< N^OMEGAl^AA 


T^N ftftOA 

SF7 = 5>.0*TF7/f 1 .0^TF7A*7l 


ISN 0096 

TG2=CN^2.0*DME6A >*A4 


—ISN OJQSL6 

SG2s2.Q«TG2yi 1 .04TG2442> 


ISN 0097 

TH=2.0*CIMEGAAA4 


-ISN JiOQfi 

SH=2.04TM/M .OtTH**21 


ISN 0099 

TFl -2 .04< N-OMEGA l*A4 


tS»J ftlftft 

SF1=2.0*TF1/I 1 .0^TF1442) 


ISN Old 

TG»N4AA 


-XSN 0102 - 

SG=2 . 0* TG/( 1 . 0+TG442 > 


ISN 0103 

TGI N*2.0*0MEGA )♦ A4 


.-LSN_01fl.A 

SGlF2.0*TGl/tl.0^TGl442) 



C.. ...COMPUTE THE TERMS. 


XSR-jQiaS T IfcQ .5aQ»PJ^jat*SEI^ 


ISN 

0106 

T2=2.0AP**4FK4*4*SF 


ISN 

0107 

T3==0.5004KF*8*SF2 


ISN 

OIOS 

T4=P*464K442*SGI 


ISN 

010 9 

TS=P4424^K4*24( I P442-K442 ) ♦♦2) 4SG 


ISN 

01! 0 

T6=P4*2FK»*6AS62 


LSN 

01 1 1. 

T7*P4484SF1 


TSN 

0112 

T85K4484SF2 


ISN 

011-3 - 

T9s4.0*P4*6*K*424SGl 


ISN 

0114 

T10=4.04P*A24K*464SG2 


LSN. 

Oils 

IIJ^OjtO*P444*K*«*SH . 


ISN 

0116 

T12=0.5D04P4474K!ltSFl 


ISN 

jQ117_ . 

tl3*P4434K*A34SP 


ISN 

0118 

T14=0.500 4PFK4474SF2 


1 SN 

.0.119 

T1 5= 1 .500*P**3*K*^34< P442-K442I ASH 


ISN 

0120 

T16=0.500AP*A5iK4( PAA2-3 . 0AKAA2 >ASG1 


ISN d2 1 

Tir=0.5D0*P*K4UP*42-KAA2)AA2l*5G 


ISN 

0122 

Tie=s0.5D0AP#'KAA5*< 3.0APAA2-KAA2 |ASG2 


ISN 

0123 . 

T19^0.5DOAP*A7AKASFl 


ISN 

0124 

T20*P4*3AK443*f PAA2-KA42 JASF 


_1SN_ 

012S 

T21=0.500AP*K*A7*SF2 


ISN 

0126 

T22=0.S00*PA*5AK *1 PAA243 .0AK4A2 )4SG1 


ISN 0127 

T 23=0 « 5D0AP* KA | C P* A2-KAA 2 | AA3 > A5G 


ISN 

0I2B 

T24=0.500AP*KAA5A< 3. OAP» A2AK4A2 1 ASG2 


ISN 

3 

T2S=1 .5D0APAA3AKA43ASH 


C44444FH0 SECTION. 
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ISN 0177 PJ-OJ/J 

ISN 0178 TF <RJ - 0*0) 49f48f4fl 

ISN 0 - 49 R J=-t9 J ^ 

ISN Oiec 49 CONTINUE 

ISN 0191 IF tRl - CRIT) 24*^4^35 

ISN 0182 34 IF <RJ - CRIT) 19»19«3S 

ISN Q1B3 35 OEL^IsDEJjry^O 

ISN 0184 IF (NA •EQ* IJ GO TO 36 

-JSN-Qlfi6- UL_JJsLA._.£aA -21- GC.. XD_37_^ 


ISN 018 8 


19 CONTINUE 
C***44END SECTION. 


C***44THIS SECTION HALVES THE STEP SIZE WHEN NCI OR NC2 EQUALS 1. NCMECK 
C keeps track of MHETHER .IHg,_.STf,P._STZ^ IS HAt,VEP FOR THE CHECK # 


ISN 

01B9 


IF CNB .EQ. 2) 

GO 

TO 

41 

ISN. 

OlSLl , - 



jiQ.ra.A2 

— - 


.. 

ISN 

0192 

41 

NCHECK»-NCHECK 




ISN 

019 3 


IF (NCHECK *E0 

.11 

_S0 TQ-.4a.--. 

1 SN 

0195 


GO TO 42 




ISN 

0196 

43 

OELT^DELT/2idJ_ 




, 

ISN 

0197 


IF CNA .EQ. 1 » 

GD 

TO 

36 

7 

0 19Q 


TF IMA »EQ. 2) 

60 

TO 

37 

ISN 

0201 

42 

CONTINUE 





rSNL 

ISN 

_1SN 

ISN 

ISN 

ISN 

ISN 

ISN 


Q.2IL2L. 


..C** SE CTI ON. 

CA4#*4iTHIS SECTION INCREMENTS THE IMPORTANT QUANTITIES. 
&PSI =0.1__i — 


020 3 
JD20A- 
0 20 5 
020 6 
0207 

020a 

02C9 


I-I 4 DI 
J^J » QJ 


XII =1 80. OAI/3. 14159 
X 180.0* J/3 • 1 4 1 59 
PSI*PSI 4 01 + DJ 
XPS 1 *18 0. 0 *P S-I >^3 . 1 4 1 59 _ 
LEsl.£+r>LE 


.I.SN. ii.210 LM^LM+Ca-M „ . - - ^ 

ISN 0211 N=LE/C 

ISN 0212 0MEGA=A3/(XI**3) . 1- 

ISN 0 213 Xl=LM/e 

.1SN -021A TETtPELT.- - 

C*****END SECTION. 

C^.jt.jtNJ'sNUMBE^^ iterations DONE SO FAR IN A RUN. 

ISN 0215 NT=NT+1 

C* * *4* THI S_ SECTJt 04 COMOUTES THE CH^N^ J IN THE LIMIT OF 

C INFINITE VlFCnSlTr {DARWIN IRSo" PAGE 317,) 

C-f .jlCPMPUTATION J.S XI 0001 TO AVOIO OIVIDING ^ B^p . _ 

C.....TH1S SECTION IS USEO ONLY IN THE CONSTANT VISCOSITY PROGRAM. 


ISN 0217 
ISN 0215- . 
ISN 0219 

ISN 022-0 

ISN 0221 

- ISN 0222 

ISN 022 3 
_ .02 24 - 

ISN 0225 
. I SN^ 022i&- — 
ISN 0227 


GO T O 33 

IF <X1 - l.COOlDOl 31.32.32 

^ OIJL«P.OOO 

DJJ*0.0D0 

GO._TO .33 , 

32 L=C*OMFGA/(Lf-LMi 

_ R1=4.0*L*( 1 .0-"LI/C1 .0-2.0»L I 

R2=0.5D04U.0/CLT-LM) ♦ 1.6/LM)*U.O + Rl ) 

R 3st-0.5D0*( 1 ._p / CLT- LM) » I . 0/LMI*< 1.0- Pl > 

Dt 1=*'DLM4R2*1 

0JJ«DLM*R3*1 

33 CONTINUE 

SECTION . 

C*****THIS SECTION PRINTS OUT THE NEW VALUES FOR THE IMPORTANT 
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- ■ - C - •• - 

ISN 0228 I^RITE (6«3) T # Kt #N • Xc>SI t OM^GA iVl S.TEMP.xi I ♦ X J,Dl *D J .OPS I 

-tSM- 022-^ 3 - -EO P M AT - ( 1 X 1 0 . 4 # I D 1 3 S . 1 X^O 1 OJkS.^ 1 7 ^3. 2 1 l.X«D 10*3 1 . IX * D 1 

I 2(1X. F9.3) »3(iX.D10«4) I 

C*tA%lrE^*D.-SeCT.lON*-^ - - 

C****%THIS SECTIDNJ DECIDES WHETHER CONSTANT DEl.T OR CONSTANT DLM SHOULD 
C - . eE-..USEO.*.- - 

C.* .••STATEMENT 12 GIVES CONSTANT STEO SIZE IN TIME. 

*STATEMgNTS-33. JVNCx6a- 


ISN 0230 IF INT - NO) 12.12. 13 

i SN i3^23 1 13 PELT -DABSCO .^1^7 2089004X1 4*1 2/( A 1 *0 »SQQ11tTAAXSl, 1 )— 

ISN 0232 IF <NA .EQ. 2) GO TO AO 

JSN- 4>33A GO TO -At 

ISN 0235 AO OELT=DAaSI0.007208900*XT**12/< AI*O.SOa*(T7-Ta+T9-TIO-Tl I))! 

4-Sht-^Oi3A -AX- CONT^nUE-, . . . 

ISN 0237 GO TO 14 

ISN 0238 - 12 OELTsOEI^Tl - 

ISN 023 9 14 CONTINUE 

— C4«L#T*END SECTIHN. .. 

€••.•• SHOUUD INTEGRATION BE FORWARD OR BACKWARD IN TIME? 

ISN 024 C IF. f MTIMF IJ nFl Ts-nFt T 

C*....ts NLAST EXCEEDED? 

..ISN D24 2 1-F ( NT-IsiLASTt 4 .A.. 99 

C*..*«IS XIMAX exceeded? 

ISN- -024 3 -4 - -UF -XXl-JUMAXJ. .3.^3*99- 

C^AAAATHIS SECTION COMPUTES THE TOTAL ANGULAR MOMENTUM AT THE END OF THE 
c auhu- -XC~SEPVES -AS-A- 


ISN 0244 
ISN 024F 
ISN 024 6 
- XSN 0247 

-XSN 024 8 
ISN 024 9 
ISN 02S0 


99 cc=Dcnscp$n 

LTC~QSQRT<i-EJuA2ALM**2A2#.0ALE»JJM*i:C) 
WRITE (6.11) LT.LTC 
11 FO,RMAT (///10X*l_fiMINlTJAl^ A-NG^^ 

1 .014.8) 

_iaO- -XQNT I NUE 


.-MQM*=aX 




STOP 

Em 



APPENDIX E 


COMPUTER PROGRAM FOR SOLAR INFLUENCE 


The computer program given in this appendix is discussed in Chapter IV , 
and in the comments listed in the program itself. 
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C«««««TTa$ IS T1DE4 

THIS PROGRAM INTEGRATES THE FIRST OF DARWINS (18801 EQUATIONS 

< 250 1 r j FIND I HE ANGLE B ETWEEN T HE PLANE OF T HE NOONS O RBI T AND 

C THE PROPER PLANE FOR ANY CHOSEN VISCOSITY OF THE EARTH. 

“C □ — IS T H E ANGL E BE T WE E N T H E P LAN E OF T H E LUNAR OR B I T ANO T HE 

C PROPER PLANE. 

— C DAR W IN I 1880 I IS ; 

C ON THE SECULAR CHANGES IN THE ELEMENTS OF THE ORBIT OF A SATELLITE 

- C RE V O LV I NG ABOU T A TI OA tX-Y^ D I S T OR T EP PLAN ET 

C IN 

— C SC IEN T I FI C PAPE R S BY S I R GEO R GE H OW ARD DA RW IN# V O L# 2. PP 2 0 8-382 

C CAMBRIDGE UNIVERSITY PRESS. 1908* 

— C T H E P APE R CAW ALS O BE F O UND IN 

C PHILOSOPHICAL TRANSACTIONS OF THE ROYAL SOCIETY* VOL* 171* IBBO. 

— C PP-TTS" - -8 Q T* ^ 

C 

— C***** T H6 PR OG R A M IN TE GRA TE S T HE E QUA TIO N F R O M XI«3* 0 5 TO Xl«l *Q C BO TO 
C 3*8 EARTH RADII I* 

— C***»*I F IN T EGRA T ING FR O M LA R GE XI TO SMALL K1 * T HE IN T EGRAL I S FO UND BY 
C SUBTRACTING THE NUMBER IN THE SUM COLUMN AT SMALL XI FROM THE 

— C NUMB ER I N T H E SUM COLU M N A T LARG E X I * T H E ANGL E J A T SMALL X I I S 

C THEM T>€ EXPONENT OF THE INTEGRAL TIMES J AT LARGE XI* 

— 

C**«««TH1S SGCTION DEFINES THE MOST IMPORTANT QUANTITIES* 

C*****X I IS SaR T I E AR TB ^ NOOW OT S T ANC E / REF ER E H CE OIS T AN C E I * 

C***«*0X1 IS THE CHANGE IN XI* 

— C***.*D ZE RO I S T H E REFERE NC E D I S T ANC E t H E R E I N UN IT S OF E A RT H " R A DII . 

C ANO WHe^E Na2FOMEGA* 

— C ****. D S» SQ R T (O t eR O> * C^ZER O - S MALL A PBS * 

C.**«*N IS ThE ROTATIONAL ANGULAR VELOCITY OF THE EARTH IN 10«4>-4 <SEC 

— C (t * E * MO L Tt PL Y TflE ^ VALU E G I V E N I N T H E PRO G R A M B Y IO»F>>» TO G ET THE 

C VALUE IN C6S UNITS*) 

— C ***** Q MgSA IS THE O R B IT AL AN G ULA R V E L OCI T Y OF T H E MOO N I N UN IT S OP 

c iOFF-4 /sec* 

— C ***4«V -I S I S THE E A RT HS V I S C OS ITY I N UN I TS O F I0»»I6 C GS* 

C«****LT IS THE TOTAL ANGULAR MOMENTUM OF TH- SYSTEM IN UNITS OF 

~~C 1 0 » »40 C OS* — — — ■ — — 

C***«.C IS TIC MOMENT OF INERTIA OF THE EARTH IN UNITS OF 10«A44 CGS* 

_C C=SMALL X »B * — — 

C IS SQRTfCBlG G^SMALL A)/(0IG M4-SMALL M))«(B1G M>«(SMALL MI«DS IN 

— C U N ITS CP LM« B* XI» 

C«***«A2 IS 19/(2P(SMALL G)4(SMALL AIMSMALL W>)IN UNITS OF 10**-12 CGS* 
C***..A3 IS Sa RTtee iB QI » ( B W M» S M ALL -M )/ C S M ALL A» »» 3)/ BS*» 3 i n u n it s o f 
C lOF^-4 SEC* OMEGAsA3/(Xl*4>3)* 

-■ END - S E C T I ON * 

C 

— C- * - »-***S C7M ^ OF D AR WI NS N OT A T I O N* 

C.***.C«ZERO ::REFERENCE DISTANCE 

C.***. OME QA- "Z g R OxO M EQA A T- C^^ Z ERO 

C*»***SMALL K»CF(OMEGA-ZERO)#(C-ZERO)/(B1G G)*CBIG MI^CSMALL M) 

C 4V -y 4.TAU " Ze R0 =(3/ 2l» € B I0 B 1 »(SM ALL N >/(C - ZEW »> 3 — 

C« ••••GOTHIC SMALL G-C 2/51 «( SMALL G»/( SMALL A) 

C **»** BI G G = un i v er sal GRAV I TA TI ONAL C O NSTA N T ~ 

C***..SMALL AiRAOIUS OF THE EARTH 

-c ** ~ * - *.B i G m ^mss O F t he e a rt h 

SMALL N^MASS OF THE MOON 








1U4 


THE “EAffTHST SURFACE 


IS THE DENSITY OF THE EARTH 

C*«*««WE HAVE SUBSTITUTED BIG G FOR DARWINS MUABOVE* 

C 

C. • • « «CORRES RONDENCE BETWEEN OUR DOTATION AND DARWINS* 
C OUR XI IS DARMINS OREEK LETTER XI* 

C tP - T All RRTNE ^ “ 

C T2ER0 - TAU ZERO FOR THE MOON 

“C GOTHIC SMALL S 

C LDA -- CREEK LETTER LAMBDA 


C T - GREEK LETTER TAU 

X. PT~- GUfHIC SHALL ~M 

C Kl - KAPPA SUB-1 * K2 - 

c 


KAPPA SUB-2 


ISN 0002 DOUBLE PRECISION DXI • VIS* X! .DZEROtCtLI. DS •AZ.TP^TZERO* GG« A3* B* A4 • 

” l“^*OMEeA *^iLDA^E*T*M*Y*tG*56* rGl«sGI •IF*SF*tFi *SF1 • ALPHA* A* BETA# “ 

2 BL* ALPHP*AP«8ETAP*eP*CAM*DELTA*TERM«Kl *K2 *X1 « X2 * X3.Z1 *X4*X5*X6* 

3-Z2*OLOGJ 

ISN 0003 DXls0*C025D0 

TWSOXni DZERO= 3* B3387305D0 ~ 

ISN 0005 C«e*llD0 

TSN 0006 LT=3AV^0 

XSN 0007 OS==DSORT<OZERO» 

T5N”O0SB A^=Z.756D0 ” — ^ 

ISN 0009 TP-5*9466920-l4 

I S^TnCTID ™ rZERO* 0* 284T4D^/TDZEROTNr3r ' 

ISN 0011 GG^6»1 S&662D-7 

ISN 0012 ' A3= TZ* AOTBSO 0/IDZERD AOS I 

ISN 0013 B-3*6B1701004DS 

* • •l^VTS ~rS THE NUMBER~OF“ VI SCOSI TlES“TO"BE"l?EAD^ 

ISN 0014 READ (5*4) NVfS 

TSN 0015 A FORMAT C ISI ‘ ^ 

ISN 0016 DO 3 J«tfNVfS 

C* • *‘i *VTS^ S THE" C HD SEN VI SCOSTTY OF THE “EAPTH* 

C**««*READ IN THE VISCOSITY* 


TSN oorr “REA0-r5*“5T“VrS 

ISN OOia 5 FORMAT (010*5) 

TSN OOrV WRITE (6*1) ’ — ‘ ^ “ “ ‘ “ 

ISN 0020 1 FORMAT ( IHl ) 

ISN “0021 write I6T61 YIS 

ISN 0022 6 FORMAT (///* 29X* 1 OMVI SCOSITY« *010*5 *1 X« lOHl 0**16 CGS*///) 

ISWVOZT WRirr-TSTTI " 

15N 0024 7 FORMAT ( ///* |SX*2HX1 • 1 1X*SHDL0GJ »12X«34SUM«/// ) 

ISN 0025 AA=VIS*A2 

C«*«**EVALUATE THE INTEGRAL* 

^5N 0026 S»0Vt) — _____ 

XSN 0027 XI»3*9€D0 

— rsiT 002B ™ocr T r»i‘* 1 tee 

ISN 0029 GMEGAV A3/(XI**3) 

ORBITAL ANGULAR MOMENTUM ^ 8«XI* 

CTTi i i *N^S trtWPUTEtr BT”T1SSU«TNC THE ROON STATS IN THE eoUATORlAL PLANED 

ISN 003D N*(LT-B7XI )/C 

ISN 0031 CUAaOMEGA^N — 

ISN 0032 E«(0«5O0)*CN**2)*l*0D-B/(GG) 

tSN 00 33 r « r ZERO/ I XI * *6) 

ISN 0034 M«C*N/(B*XI} 
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^SN~<m35 


ISH 0037 
r s N DOJa 
ISM 0039 
T5N' 0040 
iS^ 0041 


ISN 0043 


ISN 004S 

I5N 004 or 


-ycTP / T 

•HERE the tangents AND SINES OF THE LAG ANGLES ARE COMPUTED# 
rG»NTA4 ' 

SGs2«0«TG/Cl *0 * TG442) 


SG132. 04TG1/4 1«0 ^ TG1442I 


rSN 0047 
rSN— 00*8 — 
ISN 0049 
TSM 0 050 — 
ISN 0051 
T 5N 00 52 — 
ISN 0053 
ISN 0G54 — 
ISN 0055 



TPlB^^O «N«A4 

SF=2#04TF/f 1 *0 ♦ TF**2) 

TF 1 9 2 . OV tTT^^UNEG AT 

SF1 = 2.04T«^1/<1*0 ♦ TF1442) 

THERE TIC I EH MS IN" THF ETOATI ON ARE XimPPTED# 

STR A I GHT FOR WARD* 

» Y4(1^07 r 2 •G*LDA4en 

ASM 

BETA^l *0 4 V — 

BL^l *0 

ALPMP ^ W * f r * C 3Ttr/TgTTr» L D A9 ET l-T 2* 0 » C I*0» "rr9zr~ 

C2*0«C 1*0'I'Y442> 4 7«04M) 

8PF-Il^ + Y*42 4 6«04M| 

GAWa< 0 *5D0 )» W * t S F A-SGI 4SGI/ SF 1 

DELTA:=r(SFt4SGl-SG-2*0*Y4SG4<Y442>4$F)/(2*0«SFt I 

TERM»( 0 . 5 0 0 * W» » C 2* 0 4 C I * 0 4 Y) 4 r S G >2* 04SGX ) /S F l 

Kl*C-ALPHA-BETA-OSQRTf < ALPHA-BETA) 4*244 *0*A48L I) /2«0 

X1*-IK 14AL PHA)4( ALPHP-BETAP > 


Ttr*Mrr 


ISN 0061 


ISN 0063 
I SN 0064 


X3«-BP4A 
2ts1 

X4SGAM *< K24ALPHA ) 

XSg B EL T A»CK14AL P MAI 

X6-TERM 


^-SN"T)tI6S— 

-rSN 006S- 
ISN. 0067 
“TSN 0068 
ISN 0069 
“tSN «70“ 
ISN 0071 
-fSN 0072 - 


C**»*«DLOGJ IS THE CHANGE IN LOG J AT ANY ONE STEP* 

D L t T G Js <2t 42 2 )»B* DXt / ^ (C»N l — 

C*****SUM IS THE SUMMATION OF THE DLOGJ S fI«E* THE INTEGRAL.) 
S-3^ » O L OG^J ^ ■■■■■ 


WRITE <6«6I X1«DL0GJ*$ 


XX=XI-OXI 
" CONT INUE 
STOP 

T!Na 



APPENDIX F 


ERRATA FOR GOLDREICH (1966) 

The following are corrections of misprints in "History of the Lunar Orbit" 
by Peter Goldreich as the article appears in Reviews of Geophysics , vol. 4, 
pgs. 411-439, 1966. I do not claim to have caught all the misprints; some of 
the corrections may result from my own misunderstanding; but this list should 
be of use to readers of this classic paper. 

pg. 416: Equation (9) should read: 

' ‘ cx>s S = cos^ ~2 cos (4*' ~ u) + sin^ y cos (3>' + u) ’ ’ 

pg. 417: Equation (12) is derived from Equation (7) by using the approximatior 
cos I = 1 “ ;flV2 . 

If this approximation is not used, then the e3q>ression in braces in 
Equation (12) will read 

‘ cos 2 1' cos 2 u 

3 3 3 

+ cos I (sin 2$' sin 2u) cos 2<I>' + -g /S^ cos 2u ’ ’ 

The derivation of Equation (13) is still permissible, since the terms 
containing (cos 2$' cos 2u) and (sin 2$' sin 2u) are periodic, so long 
as ^ u. 
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pg. 417: Equation (14) and the line below it: How M/)H enters into the discus- 
sion is not apparent (to me). 


pg. 419: Equations (19) should read: 

- ilMi) - .. 

dt ~ dt 


f * 


Equations (22), (23), (43a), and (43b) likewise need parentheses around 
the whole quantity appearing in the differentiation operator. 

pg. 419: Read GM " for " fj. " in the equation for Kj in Equations (21). 

pg. 420: Two lines above Equation (29): 

” ... be derived by multiplying 2 (a ■ c) -jr into equation 25, 

2 (b • c) into equation 26, 2 (a • b)L into equation 27 ... ” 

1 JR 

pg. 423: Equation (41): Here " ^ ^ " should be substituted for ” ^ ” on the 

right side. Clearly the author ignores the " m '• since ^ « 1. 


pg. 423: The fifth line down from Equation (40) should read: 
" of (39) having ..." 


pg. 425: The next to the last line should read: 

K Q 

" ... Next, dotting jj into equation 43a and into ... 

pg. 426: The line above Equations (51) should read: 

" Using (1), (2), (6), and (21), we observe ... " 
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pg. 426: The first equation of (51) apparently uses 

dKi _ 1 da® _ ^ 1 d(ca®) _ 2Ki 

dt dt Ca^ dt H dt 

dC ^ 

This implies « C 


pg. 426: The equation above (55) resolves vector T along two independent sets of 
orthogonal coordinates; this procedure is ambiguous. Equations (69) - 
(74) show that the expression for T is really the sum of the lunar and 
solar torques, with the/first three vectors being the lunar torque re- 
solved along (6 j, e^, 63 ) and the last three vectors being the solar torque 
resolved along (f fj , f 3 ) • 

z ^ a 

pg. 427: Equation (58): MacDonald (1964) has sign q' = sign using 

(1 - z2)i/2 

Goldreich's notation. 

pg. 428: The right sides of the last two equations of (63) should read: 

. . 2 mA 

“ . . . ^ q B (q) sin 2 S ” 

TT a” 

“ = -^^q'F(q) sin 2S ” 

TT ar 

Confusion arises here because Goldreich corrects errors in Equa- 
tions (42) and (44) of MacDonald (1964), but inadvertently includes "n" 
in the last two equations of (63). I must confess that I do not know if the signs 
of the two equations in my correction are right, since they depend on 
MacDonald^s derivation, which I could not follow in places. 
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pg. 428: The author uses slightly different notation from Kaula (1964) in Equa- 
tion (65); " m* ” is brought outside of " B„ " and written e3q>licitly in 
Equation (64). 

pg. 428: Following the notation of Kaula (1964), ” q ” has been set equal to zero 
in Equation (66). 

pg. 429: The right sides of Equations (67) and (68) should both be multiplied by 
" m " (the lunar torque) or ")R " (the solar torque). 

pg. 429: The second term of Equation (69) and of (71) should each be multiplied 
by " kj 
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TABLE 1 


The angular speeds, phase lag angles, and amplitude factors for the seven 
tides are given. Adopted from Darwin (1880). 
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TABLE 1 


Angular 

speed 

2 (n - fi) 

2 n 

2 <n + n) 

n - 2n 

n 

n + 2Q 

2 n 

Phase lag 
angle 

2fi 

2 f 


Si 

g 

^2 

2 h 

Amplitude 

factor 

Fl 

F 

Fa 

Gi 

G 

Gj 

H 





112 


TABLE 2 

The critical angle for various viscosities is given, 
from Cq, where sin 2gj is zero, to where sin 2g^ = ±1. 


e is the distance 
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TABLE 2 


V iscosity 
(poises) 

-/'c 

(degrees) 

e 

^0 

10^7 

8,5 

8 X 10-3 

00 

O 

2.7 

8 X lO-^' 

10'® 

0.85 

8 X 10-3 

O 

o 

0.27 

8 X 10-® 

1Q2' 

0.085 

8 X 10-’ 



TABLE 3 


Summary of computer data for the curves shown in Figures 15, 16, and 17. 
The computer program itself is given in Appendix C. The column labelled 'V’ 
refers to its value when the moon is at or near 3.83 earth radii distance from 
the earth. The column labelled '’Time" gives the time required for the moon to 
move from 3.83 earth radii to 10 earth radii. The quantities At, NQ, , and 
CRIT are explained in Chapter m, Section C. 





TABLE 3 


V iscosity 
(poises) 


0 

(degrees) 


Time 

(Years) 


At 

(sec X 10'®) 



10^® 

3 

570 

2.5 X 10-® 

600 

5 X 10-“* 

10'" 

3 

3600 

2.5 X 10'® 

600 

5 X 10-'* 

lo” 

3 

3,6 X 10^ 

2.5 X 10*'* 

600 

5 X 10-'* 

lo'® 

2.68 

3.6 X 10® 

5 X 10“'* 

600 

2.5 X 10-'* 

lo'^ 

0,85 

3.6 X 10® 

5 X 10-3 

600 

2.5 X 10-“ 

10^° 

0.268 

3.6 X 10^ 

5 X 10'^ 

600 

2,5 X 10-“* 

10^' 

0.085 

3.6 X 10® 

2.5 X 10-2 

2400 

2.5 X 10-'* 

10^8 

1 

3.6 X 10® 

5 X 10-'* 

1200 

1.25 X 10-'* 

10'» 

2 

3.6 X 10® 

5 X lO-** 

1200 

1.25 X 10-'* 

lo'® 

3 

3.5 X 10® 

5 X 10*'* 

1200 

1.25 X 10-'* 

lO^i 

1 

3.6 X 10® 

1 X 10-2 

1200 

1.25 X 10-'* 

10^' 

2 

3.5 X 10® 

1 X 10-2 

1200 

1.25 X 10-'* 

10^' 

3 

3.5 X 10® 

1 X 10-2 

1200 

1.25 X 10-'* 



TABLE 4 


The important quantities used in this work are listed. 
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TABLE 4 


Symbol Description 

a mean radius of the earth 
Lm/^ 


b 

c 

Co 

e 

2fi 

g 

g 

3 

Si 

i 


m 

n 

r 

t 

X 

C 


earth-moon distance 

3.83 earth radii 

eccentricity of the lunar orbit 

lag angle of the tide with speed 2 (n - D) 

gravitational acceleration at the surface 
of the earth 

lag angle of the tide with speed n 

2 g. 

5 a 

lag angle of the tide with speed n - 2D 

angle between the invariable plane and the 
earth's equatorial plane 

angle between the invariable plane and the 
moon's orbital plane 


^ ^0 Cq 


G Mm 

mass of the moon 

angular velocity of the earth 

radial distance measured from the center 
of the earth 

time 

c - Co 

polar moment of inertia of the earth 


Numerical value 
6.37 X 10* cm 
7.21 X 10'»* 

sec 

2.44 X 10® cm 

980.7 cm/sec^ 

6.16 X 10“’ sec~’ 


1.14 X lO"* sec 
7.35 X 10^* g 


8.11 X 10^^ g-cm^ 
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TABLE 4 (Continued) 


Symbol 


Description 


Numerical value 


G universal gravitational constant 

I angle between the ecliptic and the earth's 
proper plane 

angle between the earth's proper plane 
and equatorial plane 

J angle between the moon's proper plane 

and orbital plane 

angle between the ecliptic and the moon's 
proper plane 

Lg rotational angular momentum of the earth 

L„ orbital angular momentum of the 

“ earth-moon system 

Lt total angular momentum of the earth- 
moon system 

M mass of the earth 


R disturbing function 


T 


€ 


absolute temperature of the earth 



2 Cq C 




K 


19 V 
2 g a p 

Sin 1“ ( i + j) 


X il/n 





6.67 X 10-® 


cm® 

g-sec2 


34.2 X 10^® 

sec 

5.98 X lO” g 


,.1 '"o 

7.98 X 10^^ — 

V 

2.76 X 10'^® V 
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Symbol 

7T 

P 

a 

T 

V 

-A 

n 


TABLE 4 (Continued) 

Description Numerical value 

cos I (i + j) - 

density of the earth 5.5 g/cm^ 

displacement of the earth’s surface — 

^ ^ — T. = 3.37 X 10-1° 

c3 ^6 

viscosity of the earth — 

angle between the moon’s orbital plane ~ 

and earth’s equatorial plane = i + j 

critical angle = j/sin 4 f ^ — 

orbital angular velocity of the moon — 


sec~^ 



FIGURE 1 


The earth and its attendant tidal bulge is shown on the left and the moon on 
the right in the figure. The moon orbits in the equatorial plane of the earth 
in the same direction that the earth rotates. No friction is present. 

Friction is present. The diagrams are not to scale. 




(a) 






FIGURE 2 


E is a vector normal to the ecliptic, P normal to the proper plane of the 
satellite, and M normal to the plane of the satellite's orbit. M sweeps out a 
cone about P . J is the angle between M and P, and is the angle between 
E and P. 

I is normal to the invariable plane of the planet-satellite system, A normal 
to the planet’s equatorial plane, and M normal to the satellite's orbital 
plane. A and M sweep out cones about I when solar influence is negligible, 
with all three vectors lying in a single plane. 
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FIGURE 3 

Figure 7 of Goldreich (1966), showing the inclination of the moon's orbital 
plane to the ecliptic. Precession of the lunar orbit causes the inclination to 
oscillate between the two branches of the curve. 
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FIGURE 4 

— * — ♦ 

(a) E is normal to the ecliptic, P normal to the moon's proper plane. is 
small compared to J so that the two vectors are nearly parallel and the 
inclination of the moon's orbital plane to the ecliptic is nearly constant. 

(b) J/ becomes appreciable so that the orbital plane clearly does not maintain 
a constant inclination to the ecliptic. 

(c) E lies in the surface of the cone swept out by the vector normal to the 
lunar orbit and = J. 

(d) E falls outside the cone. 


The diagrams are schematic only. 
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FIGURE 5 

The upper diagram shows the moon orbiting about the earth. L„ is the 
orbital angular momentum of the system and is perpendicular to the moon's 
orbital plane. Le is the rotational angular momentum of the earth and lies 
along the earth's axis, perpendicular to the equatorial plane. \p is the angle 
between the orbital and equatorial planes. The lower diagram shows the angular 
momentum triangle. L.j. is the total angular momentum of the system. The 
magnitudes of and are denoted by Lg and L^,, respectively. 




FIGURE 6 


The rotational angular velocity of the earth n and the orbital angular velocity 
of the rhoon Cl are shown as a function of earth -moon distance. The orbit of the 
moon lies in the equatorial plane of the earth. The dashed line is the Roche 
limit and the dotted line is the distance Cq where n = (3,83 earth radii). 







FIGURE 7 


The angular speeds of the three principal tides are shown as a function of 
earth-moon distance. Speeds n, 2 {n - n), and n - 2fi correspond to the K^, 
Mj, and O tides, respectively. The orbit of the moon lies in the equatorial 
plane of the earth. The dashed line is the Roche limit and the dotted line is the 
distance Cq where n = 2Q (3.83 earth radii). 



133 


I 

Angular 
Speed 
{sec-1 x1 O'*) 




134 


FIGURE 8 

\ - fi / n as a function of earth -moon distance. The orbit of the moon lies 
in the equatorial plane of the earth. 
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FIGURE 9 


Expression (III- 14) divided by n ^ as a function of earth-moon distance in 
the limit of low viscosity. The moon's orbit lies in the equatorial plane of 
the earth. 

Expression (III-14) multiplied by n ^ as a function of earth-moon distance in 
the limit of high viscosity. The moon’s orbit lies in the equatorial plane of 
the earth. The function is discontinuous at Cg. 
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FIGURE 10 


Sin 2gi as a function of x for large viscosities (» 10^® poises). The 
function reaches its extreme values at - e and + e . 





FIGURE 11 


The inclination ^ as a. function of x for two different initial values of \p for 
a viscosity of 10^° poises. In both cases 4’ must fall below at x = -e 
(marked by the dot with the arrow) before the moon can pass to the outer regions. 
The solid line is discussed in the text. The dotted line shows different initial 
starting conditions. The lines are not displaced for clarity. 
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FIGURE 12 


i/;2 sin 2gj, sin 4f^, and sin 2g^ as functions of x for the case of the solid 
line shown in the previous figure. Sin 2g^ Is not to scale; it is reduced by a 
factor of 10® compared to the other two functions. Sin 4fj is nearly constant. 






FIGURE 13 


dt 


for the case of the solid line shown in Figure 11. 
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FIGURE 14 

The inclination i/' as a function of x for 10^® poises for a large initial 
value of \p. The moon moves toward the earth until it reaches point D. There- 
after it moves away from the earth. must drop below the critical angle 
(marked by the dot with the arrow) before the moon can pass into the outer 
regions . 
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FIGURE 15 


The inclination \p a.s a, function of earth-moon distance for viscosities of 
10^5 , 10 and 10^’ poises for an initial perturbation of 3® at Cg (3.83 earth 
radii). 
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FIGURE 16 


The inclination i/r as a function of earth-moon distance for viscosities of 
10^® , 10^® , 10^*^, and 10*^ poises. In each case </r = at c^ - €. 
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FIGURE 17 


The inclination i/< as a function of earth-moon distance for 10^® poises 
(solid lines) and 10^^ poises (dashed lines) for perturbations of 1“, 2°, and 3® at 
Cq (3.83 earth radii). 
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FIGURE 18 

The angle i as a function of earth-moon distance for 10^® poises (solid lines) 
and 10^^ poises (dashed lines) for perturbations in of 1“, 2°, and 3° at c,j 
(3.83 earth radii). 



155 



Earth-Moon Distance 
(earth radii) 



156 


FIGURE 19 

The angle j as a function of earth-moon distance for 10^® poises (solid 
lines) and 10^*^ poises (dashed lines) for perturbations in lA of 1°, 2°, and 3° at 
C() (3.83 earth radii). 
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FIGURE 20 


The inclination. J of the moon's orbital plajie to its proper plane for various 
formulations of tidal friction. The dashed line is derived from Goldrelch (1966), 
where the three principal lag angles are equal to each other. The dotted line is 
Darwin's result for low viscosities ( «10^® poises). The upper solid line shows 
J for a perturbation in \p of 3“ at c^, (3.83 earth radii) for a viscosity of 10 
poises. The lower solid line shows J for a perturbation of 2.5° in 0 at c^ (3.83 
earth radii) for a viscosity of 10^® poises. The dashed line, dotted line, and 
lower solid line all give the present value of J at the present distance of 60 earth 
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FIGURE 21 

The point O is the center of mass of the earth and Q the center of mass of 
the moon. The earth and moon circle P, the center of mass of the earth-moon 
system, with angular velocity H . The earth rotates about the z* axis with 
angular velocity n. Vectors n and h are displaced for clarity, r and r* are the 
position vectors of the moon and mass element, respectively. 0 is the angle 
between r and r*. 
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FIGURE 22 


The position vector of the exterior point E is A. S 
a mass element in the earth. The angle between A and 


is the position vector of 
S is 'P. 
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